clear all, close all

%%______ importing data data______________

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

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

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

%calculate shear velocity
window = 5; %window for smoothing the velocity data, 5 is minimum
i=window+1;

v_shear = zeros(arraylength);

while i<length(time)-window
    v_shear(i) = (d_shear(i+window)-d_shear(i-window)) ./ (time(i+window)-time(i-window))*1000;
    i=i+1;
end

%% Plot 1 -- main

threshold = find(time == 3855);
figure(1);
subplot(2,1,1);
yyaxis left; 
plot(time,tau,'LineWidth',1.5, 'color', [0 0.4470 0.7410]); 
ax = gca;
ax.YColor = [0 0.4470 0.7410];
xlim([-inf,inf]);
ylabel('shear stress [MPa]');
yyaxis right; 
plot(time,Pp,'LineWidth',1.5, 'color',[0.6350 0.0780 0.1840]); 
ylabel('pore pressure [MPa]');
ax = gca;
ax.YColor = [0.6350 0.0780 0.1840];
ax.FontSize = 12;
set(gca, 'FontName', 'Arial');
grid on;
ax.LineWidth = 1;
xlim([3707, 11100]);
title('JY-4-PS6');

figure(1);
subplot(2,1,2);
plot(time, d_shear-d_shear(threshold(1)),'LineWidth',1.5, 'color',[0.11100 0.3250 0.0980]); 
hold on
xlabel('Time (s)');
xlabel('Displacement [{\mu}m]');
ax = gca;
ax.YColor = [0.11100 0.3250 0.0980];
ax.FontSize = 12;
set(gca, 'FontName', 'Arial');
grid on;
ax.LineWidth = 1;
xlim([3707, 11100]);
ylim([0, inf]);

set(gcf,'Position',[100 100 500 500]);

%% Tangential stiffness - unloading
% pick the data points at start and end of each stress loading
cof = 0.4; % coefficient of friction 
timepoint_start = [5520 5714 5893 6238 6520 6863 7168 7475 7781 8071 8355 8646 8922 9220 9589 9877 10190 10530];
timepoint_end = timepoint_start + 25;
dPp_start = [];
dPp_end = [];
ddisp_start = [];
ddisp_end = [];
dsigma_start = [];
dsigma_end = [];

for i = 1:length(timepoint_start)
    
    tmp1 = find(time == timepoint_start(i));
    tmp2 = find(time == timepoint_end(i));
    
    dPp_start(i) = Pp(tmp1(end));
    dPp_end(i) = Pp(tmp2(end));
    
    ddisp_start(i) = d_shear(tmp1(end));
    ddisp_end(i) = d_shear(tmp2(end));
    
    dsigma_start(i) = sigma(tmp1(end));
    dsigma_end(i) = sigma(tmp2(end));
end

dPp_diff = dPp_end - dPp_start;
ddisp_diff = ddisp_end - ddisp_start;
dsigma_diff = dsigma_end - dsigma_start;

figure(2)
plot(cumsum(ddisp_diff)*1000,abs(cumsum(dsigma_diff))*cof,'.','MarkerSize',20);
xlabel('Displacement [{\mu}m]');
ylabel({'Fault strength drop', '[MPa]'});
ax = gca;
ax.FontSize = 11;
set(gca, 'FontName', 'Arial');
ax.LineWidth = 1;
xlim([0 inf]);
set(gcf,'Position',[100 100 550 250]);
text(50, 0.75,'CSS2', 'Color','k','fontweight', 'bold', 'fontsize', 11);
%% Tangential stiffness - loading

threshold4 = find(time == 3708);
threshold5 = find(time == 3760);

figure(5);
subplot(2,1,1);
title('JY-4-PS6');
subplot(2,1,1);
yyaxis left; 
plot(time,tau,'-','LineWidth',1.5); 
xlabel('Time (s)'); 
ylabel('Shear stress [MPa]');
yyaxis right; 
plot(time, d_shear*1e3,'LineWidth',1.5); 
ylabel('Displacement [{\mu}m]');
ax = gca;
ax.FontSize = 12;
set(gca, 'FontName', 'Arial');
grid on;
ax.LineWidth = 1;
xlim([time(threshold4(1)), time(threshold5(end))]);

subplot(2,1,2);
plot(d_shear(threshold4(1):threshold5(end))*1e3,tau(threshold4(1):threshold5(end)),'k','LineWidth',1.5); 
hold on
xlabel('Displacement [{\mu}m]');
ylabel('Shear stress [MPa]');
ax = gca;
set(gca, 'FontName', 'Arial');
ax.FontSize = 12;
grid on;
ax.LineWidth = 1;
xlim([d_shear(threshold4(1)), d_shear(threshold5(end))]*1e3);


set(gcf,'Position',[100 100 500 500]);