%% TEMCO 0.25 mL/min (60 grit-roughed) *Updated 1/8/24

%Matlab script for importing and plotting data from Temco cell, calculates shear stress and normal stress for a 60 degree
%fault inclination. 

% !!! axial displacement must run from 0 upwards !!!
% otherwise piston friction is not calculated correctly

clear all


%%______ importing raw data______________


raw = readmatrix('/Users/mattroseboom/Documents/Penn State EME/Research/TEMCO/Injection Moment Experiments/MR1_Granitoid_CR0p25_1.txt');


%dt = raw(:, 1) * 1e-3; % sec
data.PpB = raw((1:end), 4) * 1e-3; % MPa pump B (pore pressure)
data.PpA = raw((1:end), 3) * 1e-3; % MPa pump A (confining)
data.PpC = raw((1:end), 5) * 1e-3; % MPa pump C (axial)
data.VpB = raw((1:end), 7); % volume remaining in B
data.VpA = raw((1:end), 6); % volume remaining in A
data.VpC = raw((1:end), 8); % volume remaining in C
data.QpB = raw((1:end), 10); % mL/min
data.QpA = raw((1:end), 9); % mL/min
data.QpC = raw((1:end), 11); % mL/min
data.d_axial = raw((1:end), 16); % mm
data.extP = raw((1:end), 17); % Psi

data.PpB = data.PpB - data.PpB(1);
data.PpA = data.PpA - data.PpA(1);
data.PpC = data.PpC - data.PpC(1);


data.time = datetime(sprintfc('%d:%d:%f', [raw((1:end), 12:13), raw((1:end), 14) + raw((1:end), 15)]));
ref_time = data.time(1); %needs to store for extracting external time
data.time = seconds(data.time - ref_time);

%below is when needs to be offset for each day that has passed
idx_midnight = find(diff(data.time<0));

i=1;
while i <= length(idx_midnight)
data.time((idx_midnight(i)+1):end) = data.time((idx_midnight(i)+1):end) + data.time(idx_midnight(i)) -data.time((idx_midnight(i)+1));
disp('midnights passed (time-vector corrected)')
disp(i)
i = i+1;
end

%% import external pump data
ext_raw = readmatrix(['/Users/mattroseboom/Documents/Penn State EME/Research/TEMCO/Injection Moment Experiments/MR1_Granitoid_CR0p25_1.txt']);
data.PpM = ext_raw(:, 3) * 1e-3; % MPa pump A
data.VpM = ext_raw(:, 6); % volume remaining in A
data.QpM = ext_raw(:, 9); % mL/min

data.PpM = data.PpM - data.PpM(1);


data.ext_time = datetime(sprintfc('%d:%d:%f', [ext_raw(:, 12:13), ext_raw(:, 14) + ext_raw(:, 15)]));
data.ext_time = seconds(data.ext_time - ref_time);


%below is when needs to be offset for each day that has passed
idx_midnight2 = find(diff(data.ext_time<0));

i=1;
while i <= length(idx_midnight2)
data.ext_time((idx_midnight2(i)+1):end) = data.ext_time((idx_midnight2(i)+1):end) + data.ext_time(idx_midnight2(i)) -data.ext_time((idx_midnight2(i)+1));
disp('midnights passed (time-vector corrected)')
disp(i)
i = i+1;
end

%% merge PpM, VpM and QpM to other data using both time vectors:
% ____ has to be done before using data !!!!!!!!!!!!! ____________________
unique_ext=unique(data.ext_time);

data.PpM = (interp1(unique_ext, data.PpM(1:end), data.time));
data.QpM = (interp1(unique_ext, data.QpM(1:end), data.time));
data.VpM = (interp1(unique_ext, data.VpM(1:end), data.time));



%%   _____ adustable values: ______

d_zero = 19962; %datapoint at which the displacement is zeroed
d_start = 50583; %datapoint from where the experiment is started, jacket friction is calculated from here
data.d_axial = data.d_axial(d_zero) - data.d_axial; % change plus (T3) or minus (T1)

% external presure converted to MPa

data.extP = data.extP / 145.038;
data.extP = data.extP - data.extP(2);
data.diffP = data.PpB - data.extP


%% this is to correct offset in displacement

%data.d_axial = offset(data.d_axial,86731,87132)


clear raw ext_raw

%%______ calculate data______________

% set-up
data.axial = data.PpC; % which pump is connected to axial stress
data.confining = data.PpA; % which pump is connected to confining stress
data.faultinlet = data.PpB; % which pump is connected to fault inlet
data.faultoutlet = data.extP; % not zero if pump is connected to fault outlet

inclination = 30; % fault is at this much degrees from long axis of sample


data.axial_cor = data.axial;
data.axial_cor(d_start:end) = data.axial(d_start:end) - 0.318*data.axial(d_start:end);  % Jacket friction subtracted from axial stress
data.axial = data.axial_cor;

% calculations
% most formulas are directly copied from Ye and Ghassemi, 2018, JGR:solid
% earth "Injection-induced shear slip and permeability enhancement in
% granite fractures" 
% !! their formula for shear displacement is incorrect !!
%data.d_normal = data.d_axial*sind(inclination) - data.d_radial*cosd(inclination); %dilation, displacement normal to fault (formula from Y&G 2018)
% data.d_shear = data.d_axial*cosd(inclination) - data.d_radial*sind(inclination); % displacement along fault in mm (formula from Y&G 2018)

data.d_shear = data.d_axial/cosd(inclination);
data.d_normal = data.d_axial/sind(inclination)-data.d_shear;

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


%calculate shear velocity
arraylength = length(data.time);

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

data.v_shear = zeros(1,arraylength);

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

%% Calculating Cumulative Seismic Moment
A = 3.14.*(0.0255/2)^2 %m2 fracture area 
D = 0.0255 %m core diameter
L = 0.05145 %m core length
V = 3.14.*(D/2)^2.*L %m3 volume
Kt = 195 %MPa/m tangential stiffness (loading)

deltaU = (data.d_shear(57715:60521)-data.d_shear(57715)).*1e-3 %%m
deltaTau = data.tau(57715)-data.tau(557715:60521)
SigmaM = (deltaU.*Kt.*V).*(1e6)
%% Cumulative Injection Volume
deltaSigmaNe = data.sigmaNe(57715)-data.sigmaNe(57715:60521) % change in normal stress [MPa]
deltaU_N = data.d_normal(57715:60521).*1e-3 % displacement normal to fault [m]
Kn = 49794 % normal stiffness [MPa/m]
cum.vol1 = (A.*deltaSigmaNe)./Kn % cumulative injection volume [m^3]
%% Calculating Relative Permeability
%l = 50.36./(1e3) %m length of fracture
%c = l^2/(data.time(9781)) %m^2/s
%A = 981.7477./(1e3) %m^2 fracture area 
%b = 0.01./(1e3) %m aperature
%Kn = (A.*data.sigmaN(9781))./cum.vol1;
%k = ((10e-3).*l^2)./(2e9.*data.time(9781))
%E = Kn.*b;
%K = ((c.*10e-3)./E)
%% basic plots
%%
figure()
plot(data.d_shear, data.tau, 'r-')
xlabel ('Shear Displacement [mm]','FontSize', 14)
ylabel ('Shear Stress [MPa]','FontSize', 14)
%xlim ([1040 4140])
%ylim ([0 6])
legend('0.25 mL/min')
%%
figure()
plot(data.d_normal, data.sigmaN,'b-')
xlabel ('Normal Displacement [mm]','FontSize', 14)
ylabel ('Normal Stress [MPa]','FontSize', 14)
%xlim ([1040 4140])
set(gca,'FontSize',14,'ycolor','k')
legend('0.25 mL/min','FontSize', 12)
%%
figure()
plot(data.d_normal, data.sigmaNe,'b-')
xlabel ('Normal Displacement [mm]','FontSize', 14)
ylabel ('Effective Normal Stress [MPa]','FontSize', 14)
%xlim ([1040 4140])
set(gca,'FontSize',14,'ycolor','k')
legend('0.25 mL/min','FontSize', 12)
%% Friction-Shear Displacment vs Time
figure()
plot(data.time, data.muE,'b-')
hold on
xlabel ('Time [s]','FontSize', 14)
ylabel ('Friction Coefficient','FontSize', 14)
ylim ([0 1])
yyaxis right
plot(data.time, data.d_shear,'r-')
ylabel ('Shear Displacement [mm]','FontSize', 14)
legend('Friction Coefficient','Shear Displacement','FontSize', 14)
%% Slip Velocity vs Time
figure()
plot(data.time, data.v_shear,'b-')
xlabel ('Time [s]','FontSize', 14)
ylabel ('Slip Velocity [um/s]','FontSize', 14)
xlim ([8970 9145])
%ylim ([0 160])
legend('0.25 mL/min - Run 1','FontSize', 14)
%%
figure(1010)
clf
tiledlayout(3,1)
ax1 = nexttile;
    plot(data.time, data.Pf,'b-')
    ylabel ('Pore Pressure [MPa]','FontSize', 14)
    xlabel ('Time [s]','FontSize', 14)
    xlim ([8737 9200])
ax2 = nexttile;
    plot(data.time, (data.d_shear-3.03),'r-')
    ylabel ('Shear Displacement [mm]','FontSize', 14)
    xlabel ('Time [s]','FontSize', 14)
    xlim ([8737 9200])
ax3 = nexttile;
    yyaxis left
    plot(data.time(57715:60521), SigmaM,'o')
    ylabel ('Modulated Seismic Moment [N.m]')
    hold on
    yyaxis right
    plot(data.time(57715:60521), cum.vol1, 'r-')
    %ylim ([0 2e-8])
    ylabel ('Cumulative Injection Volume [m^3]')
    xlim ([8737 9200])
linkaxes([ax1,ax2,ax3],"x")
%%
figure()
plot(data.Pf, 'r-')
xlabel ('Index]','FontSize', 14)
ylabel ('Shear Stress [MPa]','FontSize', 14)
%xlim ([1040 4140])
%ylim ([0 6])
legend('0.25 mL/min')
%%
figure()
plot(data.time, data.PpC,'m-')
hold on
plot(data.time, data.PpB,'r-')
xlabel ('Time [s]','FontSize', 14)
ylabel ('Pressures [MPa]','FontSize', 14)
yyaxis right
plot(data.time, data.mu,'k-')
ylabel ('Friction coefficient','FontSize', 14)
legend('Axial Pressure [MPa]', 'Upstream [MPa]', 'Friction Coefficient','FontSize', 14)
%%
figure()
plot(data.PpC,'m-')
hold on
plot(data.PpA,'k-')
plot(data.PpB,'r-')
xlabel ('Index]','FontSize', 14)
ylabel ('Pressure [MPa]','FontSize', 14)
%xlim ([150 950])
set(gca,'FontSize',14)
legend('Axial Pressure [MPa]', 'Confining Pressure','Upstream','FontSize', 12)
%%
figure()
plot(data.PpB,'r-')
hold on
plot(data.tau,'b-')
xlabel ('Time [s]','FontSize', 14)
ylabel ('Pressures [MPa]','FontSize', 14)
%xlim ([150 950])
%ylim ([0 4])
yyaxis right
plot(data.d_shear, 'g-')
ylabel 'Shear displacement [mm]'
%ylim ([0 16])
set(gca,'FontSize',14,'ycolor','k')
legend('Upstream Pressure', 'Shear Stress', 'Shear Displacement','FontSize', 12)
%%
figure()
plot(data.time, data.PpC,'m-')
hold on
plot(data.time, data.PpA,'k-')
plot(data.time, data.tau,'r-')
plot(data.time, data.PpB,'g-')
xlabel ('Time [s]','FontSize', 14)
ylabel ('Pressures [MPa]','FontSize', 14)
xlim ([1550 11950])
ylim ([0 12])
set(gca,'FontSize',14)
legend('Axial Pressure [MPa]', 'Confining Pressure','Shear Stress','Upstream Pressure','FontSize', 12)
%%
figure()
plot(data.d_shear, data.tau,'b-')
xlabel ('Shear Displacement [mm]','FontSize', 14)
ylabel ('Shear Stress [MPa]','FontSize', 14)
%xlim ([150 950])
%ylim ([0 4])
legend('0.05 mL/min','FontSize', 12)
%%
figure()
plot(data.time, data.PpB, 'r-')
xlabel ('Time [s]','FontSize', 14)
ylabel ('Upstream [MPa]','FontSize', 14)
xlim ([8740 9160])
ylim ([0 2])
set(gca,'YTick',(0:0.5:2.5))
yyaxis right
plot(data.time, data.extP, 'b-')
ylabel ('Downstream [MPa]','FontSize', 14)
ylim ([0 2])
set(gca,'YTick',(0:0.5:2.5))
set(gca,'FontSize',14,'ycolor','k')
legend('Upstream Pressure', 'Downstream Pressure','FontSize', 14)
%%
figure()
plot(data.time, data.sigmaN,'b-')
hold on
plot(data.time, data.tau, 'r-')
xlabel ('Time [s]','FontSize', 14)
ylabel ('Pressure [MPa]','FontSize', 14)
xlim ([1040 4140])
ylim ([0 6])
yyaxis right
plot(data.time, data.v_shear, 'k-')
ylim ([0 20])
ylabel ('Slip Velocity [um/s]','FontSize', 14)
set(gca,'FontSize',14,'ycolor','k')
legend('Normal Stress', 'Shear Stress', 'Slip velocity')
%%
figure()
plot(cum.vol1,SigmaM,'o')
xlabel ('Cumulative Injection Volume[m^3]','FontSize', 16)
ylabel ('Cumulative Seismic Moment [N.m]','FontSize', 16)
%xlim ([1.0e-9 5.0e-8])
%ylim ([1e-1 1e2])
set(gca,'YScale','log','FontSize',14)
set(gca,'XScale','log','FontSize',14)
legend('0.25 mL/min','FontSize', 14)
%%
figure()
plot(data.time(59346:60521),SigmaM,'o')
xlabel ('Time [s]','FontSize', 16)
ylabel ('Cumulative Seismic Moment [N.m]','FontSize', 16)
%xlim ([1.0e-10 5.0e-8])
%ylim ([1e-1 1e2])
%set(gca,'YScale','log','FontSize',14)
%set(gca,'XScale','log','FontSize',14)
legend('0.25 mL/min','FontSize', 14)
%%
figure()
plot(data.time(59280:65000), cum.vol1, 'b')
xlabel ('Time [s]','FontSize', 14)
ylabel ('Cumulative Injection Volume [m^3]','FontSize', 14)
xlim ([8970 9150])
%ylim ([0 6])
legend('0.25 mL/min - Run 1','FontSize', 14,Location='northwest')
legend('boxoff') 
%% save

save(['/Volumes/HD/Penn State EME/TEMCO/Injection Moment Experiments/Verasonics/','Mech_Data_MR1_Granitoid_CR0p25_run1.mat'])

