clear all, close all

%%______ importing data data______________

data = importdata(['JY_system_resistance.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);
d_axial = d_axial(1) - d_axial;

d_radial = 0; % have to fill in when radial displacement is measured or estimate from ???

clear data

%%
%%______ calculate data______________

% set-up
axial = PpC; % which pump is connected to axial stress
confining = PpB; % 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<length(time)-window
    v_axial(i) = (d_axial(i+window)-d_axial(i-window)) ./ (time(i+window)-time(i-window))*1000;
    i=i+1;
end
clear i

d_shear = d_axial/cosd(inclination); % displacement along fault in mm (formula from Y&G 2018)

% Pf = 0.5*(faultinlet+faultoutlet); % pore fluid pressure, average from inlet and outlet (formula from Y&G 2018)
Pp = faultinlet; % uniform pore pressure
sigma = (confining-Pp)+(axial-confining)*sind(inclination)*sind(inclination); % normal stress (formula from Y&G 2018)
tau = (axial-confining)*sind(inclination)*cosd(inclination); % shear stress on sample in MPa, (formula from Y&G 2018)
mu = (tau./sigma); % coefficient of friction

% fracture area
A = pi*(2.54/2)*2.54*1e-4;

%%
figure(1);
plot(d_shear*1e3,axial,'k.','MarkerSize',3);
hold on
ax = gca;
ax.FontSize = 12;
set(gca, 'FontName', 'Arial');
ax.LineWidth = 1;
xlabel('Displacement [um]');
ylabel('Shear stress [MPa]');
xlim([0 4850]);

x0 = 500;
y0 = 500;

width=500;
height=250;
set(gcf,'position',[x0,y0,width,height])
ax.LineWidth = 1;
box on