%c. Considerando la variazione del calore radiogenico con la profonditą nella crosta superiore

Ao=2.5*10^-6; % Calore radiogenico (Upper Crust)
A1=0.15*10^-6;  % Calore radiogenico (Middle Crust)
A2=0.12*10^-6;  % Calore radiogenico (Lower Crust)
A3=0.025*10^-6;  % Calore radiogenico (Upper Mantle)

% D=10*10^3; % Spessore arricchito in calore radiogenico

Qr=15*10^-3; % Flusso di calore alla base della crosta 
Qo1=60*10^-3; % Flusso di calore superficiale

K=3.0;   %conduttivitą termica inziale (Upper Crust)
K1m=3.0; %conduttivitą termica inziale (Middle Crust)
K1=3.0;  %conduttivitą termica inziale (Lower Crust)
K2=4.0;  %conduttivitą termica inziale (Upper Mantle)

load DepthCrust.dat
MagT=zeros(301,3);
MagT2=zeros(301,3);

T0=0;

for j=1:3
    
D=DepthCrust(2,j)*1000; % Spessore arricchito in calore radiogenico

z1=0:500:DepthCrust(2,j)*1000; %profondita calcolo geoterma Upper Crust (m)
z1Middle=max(z1)+500:500:DepthCrust(3,j)*1000; %profondita calcolo geoterma Upper Crust (m)
z2=max(z1Middle)+500:500:DepthCrust(4,j)*1000; %profondita calcolo geoterma Lower Crust (m)
z3=max(z2)+500:500:DepthCrust(5,j)*1000; %profondita calcolo geoterma Upper Mantle (m)

% Stima la geoterma nei diversi strati litosferici

Tz3=T0+((D.^2.*Ao)./K)*(1-exp(-z1./D))+(Qr./K).*z1; 
Tz3Middle=T0+((Qo1/K1m).*z1Middle)-((A1/2*K1m).*z1Middle.^2);
Tz4=T0+((Qo1/K1).*z2)-((A2/2*K1).*z2.^2);
Tz5=T0+((Qo1/K2).*z3)-((A3/2*K2).*z3.^2);
TempCrust=[Tz3';Tz3Middle';Tz4';Tz5'];
MagT(:,j)=TempCrust;
Depth=[z1';z1Middle';z2';z3'];

% figure
plot(TempCrust,-Depth./1000)
hold on

% Stima il gradiente termico in ciascuno strato

TGradUC=(max(Tz3)-min(Tz3))./max(z1).*1000;
TGradMC=(max(Tz3Middle)-min(Tz3Middle))./(max(z1Middle)-max(z1)).*1000;
TGradLC=(max(Tz4)-min(Tz4))./(max(z2)-max(z1Middle)).*1000;
TGradUM=(max(Tz5)-min(Tz5))./(max(z3)-max(z2)).*1000;

T1=T0+TGradUC.*z1./1000;
MaxT1=max(T1);

T2=MaxT1+TGradMC.*z1Middle./1000;
    
% Riduci la differenza di temperatura tra i due strati
Diff1=min(T2)-MaxT1;
T2=T2-Diff1+5;

MaxTmiddle=max(T2);
T3=MaxTmiddle+TGradLC.*(z2-max(z1Middle))./1000;

% Riduci la differenza di temperatura tra i due strati
Diff2=min(T3)-MaxTmiddle;
T3=T3-Diff2+5;

MaxT2=max(T3);
T4=MaxT2+TGradUM.*(z3-max(z2))./1000;

MagT2(:,j)=[T1';T2';T3';T4'];

end

figure
for h=1:3
    plot(MagT2(:,h),-Depth)
end

save Geotherms.dat MagT -ASCII

save Geotherms2.dat MagT2 -ASCII
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%