clear all
close all
clc
%% parameters
qT = 54*10^(-3);    % W/m^2;
z = 46.5;             % H(crust)---km;
zz = 200;           % km;
TT = 0 ;            % ℃;
TP1 = 1499;         % ℃;
TP2 = 1591;         % ℃;
deltaz = 100;       % m, 0.1km;
Ac = 0.5*10^(-6);   % W/m^3;
kc = 2.6;           % W/m/K;
Alm = 0;            % W/m^3;
klm = 3.5;          % W/m/K;
%% equation
% TB = TT+((qT/k)*deltaz)-(A*(deltaz^2)/(2*k));
% qB = qT-A*deltaz; 

a = []; 

  for iteration = 1:(z*10)
      
    TB = TT + ((qT/kc)*deltaz) - (Ac*(deltaz^2)/(2*kc));
    
    qB = qT - Ac*deltaz;
   
    TT = TB;
    
    qT  =qB;
    
    a(iteration) = TB;
  end

b = []; 

H = 0;
  for iteration = 1:(z*10)
    
    H = H + 0.1;
    
    b(iteration) = H;
  end  
  
c = []; 

  for iteration = 1:((zz - z)*10)
      
    TB = TT + ((qT/klm)*deltaz) - (Alm*(deltaz^2)/(2*klm));
    
    qB = qT - Alm*deltaz;
   
    TT = TB;
    
    qT = qB;
    
    c(iteration) = TB;
  end

d = []; 

H = z;

  for iteration = 1:((zz - z)*10)
      
    H = H + 0.1;
    
    d(iteration) = H;
    
  end
  
e = []; 

H = 0;

for iteration = 1:(zz*10)
    
    H = H + 0.1;
    
    Tk1 = TP1 * exp(0.00003*4.57*10000/192.4*H/33);
    
    e(iteration) = Tk1;
    
  end
  
f = [];

H = 0;

for iteration = 1:(zz*10) 
    
    H = H + 0.1;
    
    f(iteration) = H;
    
end

g = []; 

H = 0;

for iteration = 1:(zz*10)
    
    H = H + 0.1;
    
    Tk2 = TP2 * exp(0.00003*4.57*10000/192.4*H/33);
    
    g(iteration) = Tk2;
    
  end
  
h = [];

H = 0;

for iteration = 1:(zz*10) 
    
    H = H + 0.1;
    
    h(iteration) = H;
    
  end

  %% result
a;
b;
c;
d;
e;
f;
g;
h;
%% figure
plot(a,b,c,d,e,f,g,h,'m');
set(gca, 'YDir','reverse'); 
set(gca,'YAxisLocation','right'); 
set(gca,'XAxisLocation','top'); 
xlabel('Temperature  (℃) ');
ylabel('Depth (km)');
grid on


  

 

