% Calcola la geoterma con un valore del calore radiogenico costante

% Define A
A1=0.1*10^-6;
% Qo=mean(y1);%valori medi della retta di interpolazione (Qo e Ao)
Qo=70*10^-3;

z=0:500:35000; %profondita calcolo geoterma (m)

T0=0; % T iniziale
K=3.0; %conduttivitą  termica inziale
  
%a. Considerando il calore radiogenico costante, la geoterma: 
Tz1= T0+((Qo/K).*z)-((A1/2*K).*z.^2);
figure
plot(Tz1,z./-1000)
hold on
% In asssenza di calore radiogenico
Tz0= T0+((Qo/K).*z);
plot(Tz0,z./-1000)

% Create xlabel
xlabel('Tz1 (C)','FontSize',20,'FontWeight','bold');

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

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

%b.Considerando la variazione del calore radiogenico con la profonditą
% approssimato da una funzione esponenziale
% Define A0, D, Qr

Ao=[1.5:3.5]*10^-6;
D=[5:5:15]*10^3;
Qr=[15:5:25]*10^-3;

qo=Qr+(Ao.*D);

% Plot the heat generation with depth

figure

for i=1:3
    for j=1:3
A=Ao(j)*exp(-z./D(i));
plot(A,-z./1000)
hold on
    end
end

% Create xlabel
xlabel('Radiogenic heat (microW/m3)','FontSize',20,'FontWeight','bold');

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

% Plot the temperature variation with depth

figure

for i=1:3
    for j=1:3
        for k=1:3
Tz2 = T0+((D(i).^2.*Ao(k))./K)*(1-exp(-z./D(i)))+(Qr(j)./K).*z; 
plot (Tz2,-z./1000)
hold on
        end
    end
end

% Create xlabel
xlabel('Tz1 (C)','FontSize',20,'FontWeight','bold');

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

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

Ao=3.5*10^-6; % Calore radiogenico (Upper Crust)
A2=0.12*10^-6;  % Calore radiogenico (Lower Crust)
A3=0.02*10^-6;  % Calore radiogenico (Upper Mantle)

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

Qr=20*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)
K1=3.0; %conduttivitą termica inziale (Lower Crust)
K2=3.8; %conduttivitą termica inziale (Upper Mantle)


z1=0:500:12000; %profondita calcolo geoterma Upper Crust (m)
z2=max(z1)+500:500:30000; %profondita calcolo geoterma Lower Crust (m)
z3=max(z2)+500:500:150000; %profondita calcolo geoterma Upper Mantle (m)

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

figure
plot(TempCrust,-Depth./1000)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% d. Stima la geoterma in base al gradiente di temperatura di ciascuno strato
% in base alla geoterma precedentemente calcolata (esercizio c)

% TGradUC=20;
% TGradLC=15;
% TGradUM=8;

% Stima il gradiente geotermico per ciascuno strato

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

% Initialize the matrix

Z1=length(z1);
Z2=length(z2);
Z3=length(z3);

MagT1=zeros(Z1,1);
MagT2=zeros(Z2,1);
MagT3=zeros(Z3,1);

% Stima la temperatura per ciascuno strato in base al gradiente geotermico calcolato

for i=1:Z1
    T1=T0+TGradUC.*z1(i)./1000;
MagT1(i)=T1;
end

MaxT1=max(MagT1);

for i=1:Z2
    T2=MaxT1+TGradLC.*(z2(i)-max(z1))./1000;
MagT2(i)=T2;
end

MaxT2=max(MagT2);

for i=1:Z3
    T3=MaxT2+TGradUM.*(z3(i)-max(z2))./1000;
MagT3(i)=T3;
end

% Costruisci la geoterma

NewTemperature=[MagT1;MagT2;MagT3];

figure

plot(NewTemperature,-Depth./1000)