clear all
clc

%Define thicknesses of individual layers (17, 1.5, 10.5, 2.5, 0.5, 7.5, 11)
a = 17;
b = 1.5;
c = 10.5;
t = [a b c]; %thickness (meters)

%Define specific scenario perameters
    %SQS
    E = [30 80 30]*10^9;            %Young's moduli (Pa)
    nu = [0.2 0.23 0.2];            %Poisson's ratio
    C = [1.5 15 1.5]*10^6;          %cohesive strength of each layer (Pa) (SomeUsefulNumbers & Zhang, 2011 Table 6)

%     %QSQ
%     E = [80 30 80]*10^9;            %Young's moduli (Pa)
%     nu = [0.23 0.2 0.23];           %Poisson's ratio
%     C = [15 1.5 15]*10^6;           %cohesive strength of each layer (Pa) (SomeUsefulNumbers & Zhang, 2011 Table 6)


%Define non-specific scenario perameters
alpha = 10*10^-6;                   %thermal expansion coefficient (1/C)
rho = 2650;                         %density (kg/m^3)
g = 10;                             %grav acceleration (m/s^2)
N = length(E);                      %number of layers
M = 100;                            %number of iterations
lm = 0.4;                           %pore fluid factor (ratio of hydrostatic to lithostatic pressure gradient)
mu = 0.75;                          %coeff of friction (SomeUsefulNumbers 0.5-0.8)
theta = 30*pi/180;                  %angle of internal friction

%Define individual layer coupling moduli (5 in Bourne, 2003)
m1i = E./(1- (nu).^2);
m2i = (nu.*E)./(1-(nu).^2); 
m4i = nu./(1-nu);

%Define thickness averaged coupling moduli (6 in Bourne, 2003)
m1 = sum(m1i.*t)./sum(t);
m2 = sum(m2i.*t)./sum(t);
m4 = sum(m4i.*t)./sum(t);


%Define compound coupling moduli (8 in Bourne, 2003)
M1 = ((m1*m1i) - (m2*m2i))./(m1.^2-m2.^2); 
M2 = ((m1*m2i) - (m2*m1i))./(m1.^2-m2.^2);
M4 = m4*((m1i+m2i)./(m1+m2))-m4i;

%Define regional stress path
h0 = 3000;             %initial depth
hf = 500;              %final depth
Tgrad = 25;            %geothermal gradient (degrees C per km)
T0 = Tgrad*h0./1000;   %starting temperature
Tf = Tgrad*hf./1000;   %final temperature



%Sibson Normal Fault Tensor
[syyb, sxxb, szzb] = Sibson_NormalFaultTensor(h0,lm,mu);        %calculate initial stress state assuming normal faulting environment
sxxb = -sxxb*10^6;syyb = -syyb*10^6; szzb = -szzb*10^6;            %initial stress state (Pa)

[syyf, sxxf, szzf] = Sibson_NormalFaultTensor(hf,lm,mu);        %calculate final stress state assuming normal faulting environment
sxxf = -sxxf*10^6;syyf = -syyf*10^6; szzf = -szzf*10^6;            %final stress state (Pa)

%Stress and Temperature Changes for Each Iteration
dsyy = linspace(0,(syyf-syyb), M+1);                            %vertical (greatest compressive) stress path
dsxx = linspace(0,(sxxf-sxxb), M+1);                            %least compressive stress path
dszz = linspace(0,(szzf-szzb), M+1);                            %intermediate compressive stress path
dT = linspace(0,Tf-T0, M+1);                                    %temperature change

%Intialize layer stresses
sxxi = zeros(M,N);
syyi = zeros(M,N);
szzi = zeros(M,N);
syyi_isothermal = zeros(M,N);

for m = 1:M                                                     %loop that iterates through each stress change increment
    for i = 1:N                                                 %loop that iterates through each layer
        sxxi(m,i) =  sxxb + M1(i)*dsxx(m) + M2(i)*dszz(m) - M4(i)*dsyy(m) - alpha*E(i)*dT(m)./(1-nu(i)); %(7 in Bourne, 2003)
        szzi(m,i) =  szzb + M1(i)*dszz(m) + M2(i)*dsxx(m) - M4(i)*dsyy(m) - alpha*E(i)*dT(m)./(1-nu(i));
        syyi(m,i) =  syyb + dsyy(m) - alpha*E(i)*dT(m)./(1-nu(i));
        syyi_isothermal(m,i) =  syyb + dsyy(m);
    end
end

depth = syyi_isothermal./(rho*g*(1-lm));                        %depth (meters)


%Plot stress paths in each individual layer
figure,
subplot 511, plot(depth, sxxi.*10.^-6) %plot(sxxi(:,1)), hold on, plot(sxxi(:,2))
legend('Schist', 'Quartzite', 'Schist'), 
ylabel('\sigma_{xx}^i (MPa)'), xlabel('Depth (m)')

subplot 512, plot(depth,szzi.*10.^-6) %plot(sxxi(:,1)), hold on, plot(sxxi(:,2))
%legend('Schist', 'Quartzite', 'Schist'), 
ylabel('\sigma_{zz}^i (MPa)'), xlabel('Depth (m)')

subplot 513, plot(depth,syyi_isothermal.*10.^-6) %plot(sxxi(:,1)), hold on, plot(sxxi(:,2))
%legend('Schist', 'Quartzite', 'Schist'), 
ylabel('\sigma_{yy}^i (MPa)'), xlabel('Depth (m)')

subplot 514, plot(depth,(syyi-sxxi).*10.^-6) %plot(sxxi(:,1)), hold on, plot(sxxi(:,2))
ylabel('\sigma_{yy}^i-\sigma_{xx}^i (MPa)'), xlabel('Depth (m)')
%legend('Schist', 'Quartzite', 'Schist')

subplot 515, plot(depth,(syyi-szzi).*10.^-6) %plot(sxxi(:,1)), hold on, plot(sxxi(:,2))
ylabel('\sigma_{yy}^i-\sigma_{zz}^i (MPa)'), xlabel('Depth (m)')
%legend('Schist', 'Quartzite', 'Schist')

%Tensile Failure Check
TFC = zeros(M,N);
for j=1:3
for n=1:M
if szzi(n,j) >= (C(j)/2)                                  %(16 in Bourne, 2003)
    TFC(n,j) = 1;
  
end
end
end

%Shear Failure Check
SFC = zeros(M,N);                                           
for j = 1:3
for n=1:M
if (szzi(n,j)-syyi(n,j))/2 >= ((C(j)*cos(theta)) - ((szzi(n,j)+syyi(n,j))*sin(theta))/2) %(20 in Bourne, 2003)
    SFC(n,j) = 2;
end    
end
end

%Failure Criteria Matrix with Depth: 1 = tensile failure criteria met, 2 = shear failure criteria met, 3 = both failure criteria met
FC = zeros(M,4);
FC(:, 1) = SFC(:,1) + TFC(:,1);
FC(:, 2) = SFC(:,2) + TFC(:,2);
FC(:, 3) = SFC(:,3) + TFC(:,3);
FC(:, 4) = depth(:,1);