function dtemp = twolayer_evolve(t,temp,Ra_c,Ra_i,beta1,beta2,eta1,eta2,E,U)
%thermal evolution for a two layer model with 
% symmetrical top and bottom boundary layers in upper layer

 dtemp = zeros(4,1);
%dtemp is dT/dt
% 1 is temp of upper layer
% 2 is temp of lower layer
% 3 is density of the upper layer
% 4 is density of the lower layer
% 5 is [U] of the upper layer
% 6 is [U] of the lower layer

%variables are:
%   t = time (dependent var.) in Gyrs
%   temp = temperature of mantle in K
%   Ra_c = critical rayleigh number
%   Ra_i = flag for using internally heated Ra for lower layer
%   beta = nussult rayleigh number scaling exponent
%   eta1 = viscosity prefactor for upper layer
%   eta2 = viscosity prefactor for lower layer
%   E = activation energy, 0 or 525000
%   U = uranium concentration used to determine radioactive 
%       heat production in bottom layer

% constants
R_s = 6370000;  % Earth radius
R_i = 4770000;  % interface radius
R_c = 3490000;  % core radius
k   = 5.6;      % thermal conductivity
T_s = 273;      % surface temperature
c_p = 1250;     % specific heat
D = R_s-R_c;

 a1 = (R_s^3-R_i^3)*c_p/3.;
 a2 = R_s^2*k*3.1558e16;  %%  adjusting time for gyrs
 a3 = (R_i^3-R_c^3)*c_p/3.;
 a4 = R_i^2*k*3.1558e16;  %%  adjusting time for gyrs
 a5 = 3*3.1558e10*R_i^2/(R_s^3-R_i^3);
 a6 = 3*3.1558e10*R_i^2/(R_i^3-R_c^3);

Q = @heatprod;   % Q is in W/kg

 dtemp(2) = - (a4/(a3*temp(4)))*(temp(2) - temp(1))/ ...
              (  tbl(Ra_c,beta2,eta2,E,temp(2),temp(1),R_i,R_c,temp(4)) ...
               + tbl(Ra_c,beta1,eta1,E,temp(1),T_s    ,R_s,R_i,temp(3)) )  ...
            - (a4/(a3*temp(4)))*entrainment(temp(2),temp(1),temp(2),temp(1),T_s,temp(4),temp(3),eta1,E,D) ...
            + 3.1558e16*feval(Q,t,U)/c_p;
%           + 3.1558e16*22.e-12/c_p;

 dtemp(1) = - (a2/(a1*temp(3)))*(temp(1) - T_s)/tbl(Ra_c,beta1,eta1,E,temp(1),T_s,R_s,R_i,temp(3)) ...
            + (a4/(a1*temp(3)))*(temp(2) - temp(1))/ ...
              (  tbl(Ra_c,beta2,eta2,E,temp(2),temp(1),R_i,R_c,temp(4)) ... 
               + tbl(Ra_c,beta1,eta1,E,temp(1),T_s    ,R_s,R_i,temp(3)) ) ...
            + (a4/(a1*temp(3)))*entrainment(temp(2),temp(1),temp(2),temp(1),T_s,temp(4),temp(3),eta1,E,D) ...
            + 3.1558e16*feval(Q,t,7)/c_p;
%           + 0.;

 dtemp(3) = a5*entrainment(temp(4),temp(3),temp(2),temp(1),T_s,temp(4),temp(3),eta1,E,D);
 dtemp(4) = -a6*entrainment(temp(4),temp(3),temp(2),temp(1),T_s,temp(4),temp(3),eta1,E,D);


%%  subfunction to calc. volumetric entrainment rate  %%
%%  --------------------------------------------------%%
function x = entrainment(x1,x2,t1,t2,t_s,rho1,rho2,eta,E,D)

B = feval(@buoynum,rho1,rho2,t1,t2);
Ra = (feval(@rayleigh,eta,E,(t1+t2)/2,t_s,D,4500))^(.33);
%Ra = (feval(@rayleigh,eta,E,t1,t_s,D,4500))^(.33);

%x = (.2*(x1-x2)/D)*B^(-2)*Ra*(1/(1+B^-1));
 x = (.035*.35*(x1-x2)/D)*B^(-2)*Ra*(1/(1+B^-1));
%%  --------------------------------------------------%%



%%  subfunction to calc. buoyancy number  %%
%%  --------------------------------------%%
function x = buoynum(rho1,rho2,t1,t2)

%  rho1,t1 are the bigger values (than rho2,t2)
delrho = abs(rho1-rho2)/rho2;
del_T = t1-t2;
temp_b = delrho/(del_T*3.e-5);
if (temp_b < .1)
    temp_b = 1.e10
%   break
end

x = temp_b;
%%  --------------------------------------%%



%%  subfunction to calc. thermal boundary layer thickness  %%
%%  -----------------------------------------------------  %%
function x = tbl(Ra_c,beta,eta,E,t1,t2,r1,r2,rho)

% r1 = radius of upper shell
% r2 = radius of lower shell

D   = r1 - r2;  % shell depth

x = D*(Ra_c/rayleigh(eta,E,t1,t2,D,rho))^beta;
%%  -----------------------------------------------------  %%



%%  subfunction to calc. rayleigh number  %%
%%  ------------------------------------  %%
function r = rayleigh(eta,E,t1,t2,D,rho)
alph = 3.e-5;   % thermal expansivity
c_p = 1250;     % specific heat
g   = 9.8;      % gravity
kap = 1.e-6;    % thermal diffusivity

a1 = (alph*g*rho*D^3)/kap;

r = a1*(t1- t2)/visc(eta,E,t1);
%%  ------------------------------------  %%
 


%%  subfunction to calc. temp. dep. viscosity  %%
%%  -----------------------------------------  %%
function v = visc(eta,E,temp)
R   = 8.31;     % gas constant
if E == 1
   E = 525000;
end
v = eta*exp(E/(R*temp));
%%  -----------------------------------------  %%



%%  subfunction to calc. prefactor for viscosity  %%
%%  --------------------------------------------  %%
function v = viscinv(eta_pd,E,temp)
R   = 8.31;     % gas constant
if E == 1
   E = 525000;
end
v = eta_pd/exp(E/(R*temp));
%%  --------------------------------------------  %%


%%  subfunction to calc. radiogenic heat prod  %%
%%  -----------------------------------------  %%
function y = heatprod(t,u)
%function for radiogenic heat production
%t is time (Gyr) 
%u is uranium concentration in parts per billion

%radioactivity constants
lambda238 = .155;
lambda235 = .985;
lambda232 = .0495;
lambdak40 = .555;
H238  = 9.37e-5;
H235  = 5.69e-4;
H232  = 2.69e-5;
Hk40  = 2.79e-5;
ratio238 = .9927;
ratio235 = .0072;
ratio232 = 4.01;
ratiok40 = 1.28;
%t0 = 1.4359e+17;
t0 = 4.5;

a238 = H238*ratio238*1.e-9;
a235 = H235*ratio235*1.e-9;
a232 = H232*ratio232*1.e-9;
ak40 = Hk40*ratiok40*1.e-9;

y = u*a238*exp(lambda238*(t0-t)) ...
  + u*a235*exp(lambda235*(t0-t)) ...
  + u*a232*exp(lambda232*(t0-t)) ...
  + u*ak40*exp(lambdak40*(t0-t));
%%  -----------------------------------------  %%


%%  subfunction to calc. present day temperature  %%
%%  --------------------------------------------  %%
function y = temp_pd(eta_pd,D,Ra_c,beta)

%constants
k   = 5.6;      % thermal conductivity
T_s = 273;      % surface temperature
rho = 4500;     % density of upper layer
c_p = 1250;     % specific heat
g   = 9.8;      % gravity
kap = 1.e-6;    % thermal diffusivity
alph = 3.e-5;   % thermal expansivity
%q_pd = 4.e13;   % present day surface heat loss
q_pd = .0784;

Ra = alph*g*rho*D^3/(kap*eta_pd);
delt = q_pd*D/k;

y = (delt*(Ra_c/Ra)^beta)^(1/(1+beta)) + T_s;
%%  --------------------------------------------  %%

