clear all, close all
%%______ importing data data______________
data = importdata('JY_19_PS16_data.txt', '\t');
% dt = data(:, 1) * 1e-3; % sec
PpB = data(:, 3) * 1e-3 + 25e-3; % MPa pump B
PpA = data(:, 2) * 1e-3; % MPa pump A
PpC = data(:, 4) * 1e-3; % MPa pump C
VpB = data(:, 6)*1e-6; % volume remaining in B
VpA = data(:, 5)*1e-6; % volume remaining in A
VpC = data(:, 7)*1e-6; % volume remaining in C
QpB = data(:, 9); % mL/min
QpA = -data(:, 8); % mL/min
QpC = data(:, 10); % mL/min
d_axial = data(:, 15); % mm
time = data(:,11)*60*60 + data(:,12)*60 + data(:,13);
time = time - time(1);
clear data
%%
% set-up
axial = PpB; % which pump is connected to axial stress
confining = PpC; % which pump is connected to confining stress
faultinlet = PpA; % which pump is connected to fault inlet
faultoutlet = 0; % not zero if pump is connected to fault outlet
inclination = 30; % fault is at this much degrees from long axis of sample
%calculate piston velocity
window = 15; %window for smoothing the velocity data, 5 is minimum
i=window+1;
arraylength = [length(time) 1];
v_axial = zeros(arraylength);
while i