clear all, close all

%%______ importing data data______________

data = importdata(['JY_normal_stiffness_test2.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 = 90; % 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

% correction for piston friction (velocity dependent) and jacket friction (displacement dependent)
mu_piston = 0.0078.*log(v_axial.*1000).*log(v_axial.*1000) - 0.0919.*log(v_axial.*1000) + 0.5840; % piston friction
axial_cor = axial - mu_piston;

axial_cor = axial_cor - 0.2442*d_axial - 1.5835; % Jacket friction subtracted from axial stress

d_normal = d_axial*sind(inclination) - d_radial*cosd(inclination); %dilation, displacement normal to fault (formula from Y&G 2018)
d_shear = d_axial*cosd(inclination) - d_radial*sind(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

%% Plot 1 -- main

figure(1);
subplot(2,1,1);
yyaxis left; 
plot(time,d_axial*1e3,'LineWidth',1.5, 'color', [0 0.4470 0.7410]); 
ax = gca;
ax.YColor = [0 0.4470 0.7410];
ylabel('aperture change [um]');
yyaxis right; 
plot(time,sigma,'LineWidth',1.5, 'color',[0.6350 0.0780 0.1840]); 
ylabel('normal stress [MPa]');
ax = gca;
ax.YColor = [0.6350 0.0780 0.1840];
ax.FontSize = 12;
set(gca, 'FontName', 'Arial');
grid on;
ax.LineWidth = 1;
% xlim([time(threshold1(1)), time(threshold2(end))]);
title({'Normal stiffness #2'});

figure(1);
subplot(2,1,2);
timepoint = [86 212 344 468 602 731 871 1020 1129 1274 1403 1527 1634 1758 1887 2016 2117 2301 2457 2608 2774 2892 3032 3156 3285 3419 3543 3661 3806 3946];
dV = [];
dsig = [];
for i = 1:length(timepoint)
    tmp = find(time == timepoint(i));
    ddisp(i) = d_axial(tmp(end));
    dsig(i) = sigma(tmp(end));
end
ddisp1 = ddisp - ddisp(1);    % aperture change data point [um]
dsig1 = dsig - dsig(1);    % sigma change data point
plot(ddisp1*1e3,dsig1,'ko','MarkerSize',5);
hold on
ax = gca;
ax.FontSize = 12;
set(gca, 'FontName', 'Arial');
grid on;
ax.LineWidth = 1;
xlabel('aperture change [um]');
ylabel('normal stress [MPa]');

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

save('JY_16_kn_90.mat','ddisp1','dsig1')