%% TEMCO: Uniform Pore Pressure Profile


clear all


%%______ importing raw data______________


raw = readmatrix('.txt');


%dt = raw(:, 1) * 1e-3; % sec
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.extP = raw((1:end), 16); % Psi

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); %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(['.txt']);
data.PpM = ext_raw(:, 2) * 1e-3; % MPa 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);

%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));

%%   _____ adjustable values: ______

d_zero = 4895; %datapoint at which the displacement is zeroed
d_start = 5649; %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

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 referenced from Ye and Ghassemi, 2018, JGR:solid

data.d_shear = data.d_axial.*cosd(inclination);
data.d_normal = data.d_axial.*sind(inclination);
d_shear = data.d_shear; %[mm]

% Average pore pressure
data.Pf = 0.5*(data.faultinlet+data.faultoutlet); % pore fluid pressure, average from inlet and outlet (formula from Y&G 2018)

% Normal and shear 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)

%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