% This program estimates the strength and the effective elastic thickness

% Definition of the input parameters for the brittle strength
f1=0.75;
f2=3;
lambda=0.36;
g=9.8;

% Definition of the input parameters for the ductile strength
StrainR=10.^-15;
T=273;
R=8.314;

% Definition of the imput parameter to calculate ductile strength of the upper crust for each temperature values

% Quartzite (wet) (T ductile=318°C)

N1=1.9;
B1=1.26e-13;
E1=172600;

% Quartzite (dry) (T ductile=383°C)

% N1=2.72;
% B1=6.03e-24;
% E1=134000;

% Granite (dry) (T ductile=505°C)

% N1=3.3;
% B1=3.16e-26;
% E1=186000;

% Definition of the imput parameter to calculate ductile strength of the middle crust for each temperature values

% Diorite (wet)     (T ductile=433°C)

N2=2.4;
B2=1.26e-16;
E2=212000;

% Diabase (dry)     (T ductile=575°C)

% N2=3.05;
% B2=6.31e-20;
% E2=276000;

% Definition of the imput parameter to calculate ductile strength of the lower crust for every temperature value

% Mafic Granulite(?) (T ductile=730°C)

% N=4.2;
% B=8.8334e-22;
% E=445000;

% Diorite (wet)     (T ductile=433°C)
N=2.4;
B=1.26e-16;
E=212000;

% Diabase (dry)     (T ductile=575°C)
% N=3.05;
% B=6.31e-20;
% E=276000;

% Definition of the imput parameter to calculate ductile strength of the mantle lithosphere for every temperature value
% Demouchy et al., 2013 (Dry Olivine T  ductile=790°C)

StrainR1=10.^-15;
Ap=1*10.^6;
Ep=450000;
sigmaP=15*10^9;

% load TempLitho2.dat

load Geotherms2.dat

% load FinalTemperatureLitho.txt
% T0=TempLitho2;

 T0=Geotherms2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load DensityRocks.dat
load DepthCrust.dat

k=length(DepthCrust(1,:));

% Define the depth of each lithospheric layer

depthS=DepthCrust(1,:);
densityS=DensityRocks(1,:);
depthUC=DepthCrust(2,:);
densityUC=DensityRocks(2,:);
depthMC=DepthCrust(3,:);
densityMC=DensityRocks(3,:);
depthLC=DepthCrust(4,:);
densityLC=DensityRocks(4,:);
depthM=DepthCrust(5,:);
densityM=DensityRocks(5,:);

% Initialize all the matrix 

% Brittle Strength
MagBrittleStrengthComp=zeros(length(T0),k);
MagBrittleStrengthTen=zeros(length(T0),k);

% Strength (Brittle+Ductile) Lithosphere
MagBrittleDuctileStrengthComp2=zeros(length(T0),k);
MagBrittleDuctileStrengthTen2=zeros(length(T0),k);

MagPercComp=zeros(k,1);
MagPercTen=zeros(k,1);

% % Strength (Brittle+Ductile) Crust
% MagBrittleDuctileStrengthCompCr=zeros(depthLC(i).*2,k);
% MagBrittleDuctileStrengthTenCr=zeros(depthLC.*2,k);

% Integrated Strength
MagIntStr1=zeros(k,1);
MagIntStr=zeros(k,1);
MagIntStr1Cr=zeros(k,1);
MagIntStrCr=zeros(k,1);

% Coupling/Decoupling Conditions
MagVer=zeros(k,1);
MagH=zeros(k,1);

for i=1:k
 
% Definition of the vector of the brittle strength values for the sedimentay layer

strengthT=zeros(depthS(i).*2,1);
strengthC=zeros(depthS(i).*2,1);

% Calculation of the brittle strength for the sediment layer (Coulomb criterion)

depth1S=[1:depthS(1,i)*2];
j=0.5.*2:depthS(i).*2;

for j1=1:length(j)
sigmaBT=f1*densityS(i).*g.*(j(j1)./2).*(1-lambda);
sigmaBC=f2*densityS(i).*g.*(j(j1)./2).*(1-lambda);
strengthT(j1,:)=sigmaBT;
strengthC(j1,:)=sigmaBC;
end

% Definition of the vector of the brittle strength values for the upper crust layer

strengthTc=zeros((depthUC(i)-depthS(i)).*2,1);
strengthCc=zeros((depthUC(i)-depthS(i)).*2,1);

% Calculation of the brittle strength for the upper crust layer (Coulomb criterion)

depth1UC=[(depthS(i)+0.5:0.5:depthUC(i))];
m=(depthS(i)+0.5).*2:(depthUC(i).*2);

for m1=1:length(m)
sigmaBTc=f1*densityUC(i).*g.*(m(m1)./2).*(1-lambda);
sigmaBCc=f2*densityUC(i).*g.*(m(m1)./2).*(1-lambda);
strengthTc(m1,:)=sigmaBTc;
strengthCc(m1,:)=sigmaBCc;
end

% Definition of the vector of the brittle strength values for the middle crust layer

strengthTmc=zeros((depthMC(i)-depthUC(i)).*2,1);
strengthCmc=zeros((depthMC(i)-depthUC(i)).*2,1);

% Calculation of the brittle strength for the middle crust layer (Coulomb criterion)

depth1MC=[(depthUC(i)+0.5:0.5:depthMC(i))];
q=(depthUC(i)+0.5).*2:(depthMC(i).*2);

 for q1=1:length(q)
sigmaBTmc=f1*densityMC(i).*g.*(q(q1)./2).*(1-lambda);
sigmaBCmc=f2*densityMC(i).*g.*(q(q1)./2).*(1-lambda);
strengthTmc(q1,:)=sigmaBTmc;
strengthCmc(q1,:)=sigmaBCmc;
end

% Definition of the vector of the brittle strength values for the lower crust layer

strengthTLc=zeros((depthLC(i)-depthMC(i)).*2,1);
strengthCLc=zeros((depthLC(i)-depthMC(i)).*2,1);

% Calculation of the brittle strength for the lower crust layer (Coulomb criterion)
  
depth1LC=[(depthMC(i)+0.5:0.5:depthLC(i))];   
p=(depthMC(i)+0.5).*2:(depthLC(i).*2);

for p1=1:length(p)
sigmaBTLc=f1*densityLC(i).*g.*(p(p1)./2).*(1-lambda);
sigmaBCLc=f2*densityLC(i).*g.*(p(p1)./2).*(1-lambda);
strengthTLc(p1,:)=sigmaBTLc;
strengthCLc(p1,:)=sigmaBCLc;
end

% Definition of the vector of the brittle strength values for the upper mantle layer

strengthTM=zeros((depthM(i)-depthLC(i)).*2,1);
strengthCM=zeros((depthM(i)-depthLC(i)).*2,1);

% Calculation of the brittle strength for the upper mantle layer (Coulomb criterion)

depth1M=[(depthLC(i)+0.5:0.5:depthM(i))];   
r=(depthLC(i)+0.5).*2:(depthM(i).*2);

for r1=1:length(r)
sigmaBTM=f1*densityM(i).*g.*(r(r1)./2).*(1-lambda);
sigmaBCM=f2*densityM(i).*g.*(r(r1)./2).*(1-lambda);
strengthTM(r1,:)=sigmaBTM;
strengthCM(r1,:)=sigmaBCM;
end

% Calculation of the brittle strength for the sedimentary, upper, middle, lower crust and upper mantle layers

strengthT1=strengthT.*1000;
strengthC1=strengthC.*1000*-1;
strengthTcc=strengthTc.*1000;
strengthCcc=strengthCc.*1000*-1;
strengthTmc1=strengthTmc.*1000;
strengthCmc1=strengthCmc.*1000*-1;
strengthTLcA=strengthTLc.*1000;
strengthCLcA=strengthCLc.*1000*-1;
strengthTMa=strengthTM.*1000;
strengthCMa=strengthCM.*1000*-1;

BrittleStrengthComp=[0;strengthC1;strengthCcc;strengthCmc1;strengthCLcA;strengthCMa];
BrittleStrengthTen=[0;strengthT1;strengthTcc;strengthTmc1;strengthTLcA;strengthTMa];

MagBrittleStrengthComp(:,i)=BrittleStrengthComp;
MagBrittleStrengthTen(:,i)=BrittleStrengthTen;

FDepth=[0:0.5:(length(MagBrittleStrengthComp(:,i))-1)/2];
% plot(MagBrittleStrengthComp(:,i),-FDepth)
% hold on
% plot(MagBrittleStrengthTen(:,i),-FDepth)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calculation of the ductile strength in the upper crust (power-law creep)

sigmaCr2=(((StrainR/B1).^(1./N1)).*exp(E1./(N1.*R.*(T0(2:length(depth1UC)+1,i)+T))));

StrT=[strengthTcc';sigmaCr2'];
StrT1=min(StrT);
StrC=[strengthCcc';-sigmaCr2'];
StrC1=min(abs(StrC));

% Calculation of the ductile strength in the middle crust (power-law creep)

sigmaCr3=(((StrainR/B2).^(1./N2)).*exp(E2./(N2.*R.*(T0(length(depth1UC)+1+1:length(depth1MC)+length(depth1UC)+1,i)+T))));
StrTmc=[strengthTmc1';sigmaCr3'];
StrT1mc=min(StrTmc);
StrCmc=[strengthCmc1';-sigmaCr3'];
StrC1mc=min(abs(StrCmc));

% Calculation of the ductile strength in the lower crust (power-law creep)

sigmaCrl=(((StrainR/B).^(1./N)).*exp(E./(N.*R.*(T0(length(depth1UC)+1+1+length(depth1MC):length(depth1LC)+length(depth1MC)+length(depth1UC)+1,i)+T))));
StrTlc=[strengthTLcA';sigmaCrl'];
StrT1lc=min(StrTlc);
StrClc=[strengthCLcA';-sigmaCrl'];
StrC1lc=min(abs(StrClc));

BrittleDuctileStrengthCompCr=[0;StrC1';StrC1mc';StrC1lc'];
BrittleDuctileStrengthTenCr=[0;StrT1';StrT1mc';StrT1lc'];

% Calculation of the ductile strength profile in the mantle lithosphere (Dorn creep law)

cn=R.*(T0(length(depth1UC)+1+1+length(depth1MC)+length(depth1LC):length(depth1M)+length(depth1LC)+length(depth1MC)+length(depth1UC)+1,i)+T);
cn2=log(StrainR1./Ap);
cn3=(-1.*(cn./Ep).*cn2);
Dorn=sigmaP.*((1-sqrt(cn3)).^2);

FminDorn=find(Dorn==min(Dorn));
Dorn(FminDorn+1:length(Dorn))=0;

StrTM=[strengthTMa';Dorn'];
StrT1M=min(StrTM);

StrCM=[strengthCMa';-Dorn'];
StrC1M=min(abs(StrCM));

BrittleDuctileStrengthComp2=[0;StrC1';StrC1mc';StrC1lc';StrC1M'].*-1;
BrittleDuctileStrengthTen2=[0;StrT1';StrT1mc';StrT1lc';StrT1M'];

MagBrittleDuctileStrengthComp2(:,i)=BrittleDuctileStrengthComp2;
MagBrittleDuctileStrengthTen2(:,i)=BrittleDuctileStrengthTen2;

figure

plot(BrittleDuctileStrengthTen2,-FDepth)
hold on
plot(BrittleDuctileStrengthComp2,-FDepth)

% Create xlabel
xlabel('Strength (Pa)','FontSize',20,'FontWeight','bold');
% Create ylabel
ylabel('Depth (km)','FontSize',20,'FontWeight','bold')

% Estimate the Integrated Strength

dz=500;
IntStr1=sum(BrittleDuctileStrengthComp2).*(dz);
IntStr=sum(BrittleDuctileStrengthTen2).*(dz);
IntStr1Cr=sum(BrittleDuctileStrengthCompCr.*(dz));
IntStrCr=sum(BrittleDuctileStrengthTenCr.*(dz));

MagIntStr1(i)=IntStr1;
MagIntStr(i)=IntStr;
MagIntStr1Cr(i)=IntStr1Cr;
MagIntStrCr(i)=IntStrCr;

% Estimate the percentage of crustal strength

MagPercComp=(MagIntStr1Cr./MagIntStr1).*100;
MagPercTen=(MagIntStrCr./MagIntStr).*100;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the Elastic thickness of the lithosphere for all conditions (when the layers are welded together or not)
    
d4=min(sigmaCr2);
d4A=min(sigmaCr3);
d04=min(sigmaCrl);

 % Estimate the mechanical thickness of each lithospheric layer

if (d4)>10.*10^6 & (d4A)>10.*10^6 & (d04)>10.*10^6
 
    h1=length(depth1UC)./2;
    h2=length(depth1MC)./2;
    h3=length(depth1LC)./2;
    TresholdValueM=find(Dorn<10.*10^6);
    h4=min(TresholdValueM./2);
    H=(h1+h2+h3+h4);
    VerC=1;
    
elseif (d4)>10.*10^6 & (d4A)<10.*10^6 & (d04)<10.*10^6
    h1=length(depth1UC)./2;
    TresholdValueMC=find(sigmaCr3<10.*10^6); 
    h2=min(TresholdValueMC./2);
    TresholdValueLC=find(sigmaCrl<10.*10^6); 
    h3=min(TresholdValueLC./2);
    TresholdValueM=find(Dorn<10.*10^6);
    h4=min(TresholdValueM./2);
    H1=(h1+h2);
    H=((H1.^3)+(h3.^3)+(h4.^3)).^(1./3);
    VerC=2;
    
elseif (d4)<10.*10^6 & (d4A)>10.*10^6 & (d04)<10.*10^6
    
    TresholdValueUC=find(sigmaCr2<10.*10^6); 
    h1=min(TresholdValueUC./2);
    h2=length(depth1MC)./2;
    TresholdValueLC=find(sigmaCrl<10.*10^6); 
    h3=min(TresholdValueLC./2);
    TresholdValueM=find(Dorn<10.*10^6);
    h4=min(TresholdValueM./2);
    H2=(h2+h3);
    H=((H2.^3)+(h1.^3)+(h4.^3)).^(1./3);
    VerC=3;
    
elseif (d4)<10.*10^6 & (d4A)<10.*10^6 & (d04)>10.*10^6
    TresholdValueUC=find(sigmaCr2<10.*10^6);
    h1=min(TresholdValueUC./2);
    TresholdValueMC=find(sigmaCr3<10.*10^6); 
    h2=min(TresholdValueMC./2);
    h3=length(depth1LC)./2;
    TresholdValueM=find(Dorn<10.*10^6);
    h4=min(TresholdValueM./2);
    H3=(h3+h4);
    H=((H3.^3)+(h1.^3)+(h2.^3)).^(1./3);
    VerC=4;
     
elseif (d4)>10.*10^6 & (d4A)>10.*10^6 & (d04)<10.*10^6
     h1=length(depth1UC)./2;
     h2=length(depth1MC)./2;
     TresholdValueLC=find(sigmaCrl<10.*10^6); 
     h3=min(TresholdValueLC./2);
     TresholdValueM=find(Dorn<10.*10^6);
     h4=min(TresholdValueM./2);
     H4=(h1+h2+h3);
     H=((H4.^3)+(h4.^3)).^(1./3);
     VerC=5;
     
elseif (d4)<10.*10^6 & (d4A)>10.*10^6 & (d04)>10.*10^6 
     TresholdValueUC=find(sigmaCr2<10.*10^6);
     h1=min(TresholdValueUC./2);
     h2=length(depth1MC)./2;
     h3=length(depth1LC)./2;
     TresholdValueM=find(Dorn<10.*10^6);
     h4=min(TresholdValueM./2);
     H5=(h2+h3+h4);
     H=((H5.^3)+(h1.^3)).^(1./3);
     VerC=6;
     
elseif (d4)>10.*10^6 & (d4A)<10.*10^6 & (d04)>10.*10^6 
     h1=length(depth1UC)./2;
     TresholdValueMC=find(sigmaCr3<10.*10^6); 
     h2=min(TresholdValueMC./2);
     h3=length(depth1LC)./2;
     TresholdValueM=find(Dorn<10.*10^6);
     h4=min(TresholdValueM./2);
     H6=(h1+h2);
     H6a=(h3+h4);
     H=((H6.^3)+(H6a.^3)).^(1./3);
     VerC=7;
     
elseif (d4)<10.*10^6 & (d4A)<10.*10^6 & (d04)<10.*10^6
     TresholdValueUC=find(sigmaCr2<10.*10^6); 
     h1=min(TresholdValueUC./2);
     TresholdValueMC=find(sigmaCr3<10.*10^6); 
     h2=min(TresholdValueMC./2);
     TresholdValueLC=find(sigmaCrl<10.*10^6); 
     h3=min(TresholdValueLC./2);
     TresholdValueM=find(Dorn<10.*10^6);
     h4=min(TresholdValueM./2);
     H=(h1.^3+h2.^3+h3.^3+h4.^3).^(1./3);
     VerC=8;
end

MagVer(i)=VerC;
MagH(i)=H;
end

% Estimate the rigidity

YoungMod=100*10^9;
PoissonCoeff=0.25;
Rigidity=(YoungMod.*(MagH.*1000).^3)./(12*(1-PoissonCoeff^2));

figure
plot(T0,-FDepth)

% Create xlabel
xlabel('Temperature (C)','FontSize',20,'FontWeight','bold');
% Create ylabel
ylabel('Depth (km)','FontSize',20,'FontWeight','bold')