%% TEMCO Teflon Experimental Data Deduction

clear all


%%______ Importing raw data______________


raw = readmatrix('Text_File_Name.txt');


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

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

time = raw(:,11)*60*60 + raw(:,12)*60 + raw(:,13);
time = time - time(1);

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

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('Text_File_Name.txt');
data.PpM = ext_raw(:, 2) * 1e-3; % Upstream pressure pump A
data.VpM = ext_raw(:, 5); % volume remaining in A
data.QpM = ext_raw(:, 8); % mL/min

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


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

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:
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 = 3555; % Datapoint at which the displacement is zeroed (point directly after 0.5 MPa pressure step to 3 MPa)
d_start = 3951; % Datapoint from where the experiment is started, jacket friction is calculated from here (point directly at start of axial loading to peak shear strength)
data.d_axial = data.d_axial(d_zero) - data.d_axial; % Change plus (T3) or minus (T1)

clear raw ext_raw

%%
%%______ Calculate data ______________

% Set-up
data.axial = data.PpC; % Axial stress [MPa]
data.confining = data.PpA; % Confining stress [MPa]
inclination = 30; % Fault angle

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

% Displacements
data.d_shear = data.d_axial*cosd(inclination); % Shear displacement [mm]
data.d_normal = data.d_axial*sind(inclination); % Normal Displacement [mm]

% Stress calculations
data.sigmaN = data.confining+(data.axial-data.confining)*sind(inclination)*sind(inclination); % Normal stress [MPa]
data.tau = (data.axial-data.confining)*sind(inclination)*cosd(inclination); % Shear stress [MPa]
mu = data.tau./data.sigmaN; % Friction coefficient

% Teflon adjustment 
TF = 5; % Teflon width [mm]
TF_w = TF.*1e-3; % Teflon width [m]
mu_R = 0.446; % Friction coefficient for Granitoid
mu_T = 0.05; % Friction coefficient for Teflon
D = 0.0255; % Core diameter [m]
A_TF = TF_w.*D; % Teflon area

% Equivalent pore pressure
deltaP_e = (sigmaN.*(1-mu_T./mu_R));
mu_e = tau./(sigmaN-deltaP_e); % Effective friction
sigmaN_e = sigmaN-deltaP_e; % 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

%% Parameters for Seismic Moment
A = 3.14.*(0.0255/2)^2; % Fracture area [m^2] 
Kt = 3946; % Tangential stiffness (loading) [MPa/m]
U_s = (data.d_shear(6045:15115)-data.d_shear(6045)).*1e-3; % Shear displacement [m]

%% Equivalent Volume Injected
Kn = 56444; % Normal stiffness [MPa/m]
deltaVinj = A_TF.*(deltaP_e./Kn); % Equivalent injected volume [m^3]

%% Revised Seismic Moment
a = 0.03531; % Length of fracture [m]
SigmaM_e = (A.*U_s.*Kt.*a).*(1e6); % Cumulative seismic moment [N.m]