%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This code solves numerically the diffusion equation
%% dT/dt=kd^2T/dz^2 for T(z,t), emulating the solutions of both
%% Lord Kelvin and John Perry. 
%% Rei Chemke, reichemke@gmail.com, 2014
%% modifications: Eli Tziperman
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;
clear all;

%% Initialize parameters:
%% ----------------------
% depth coordinate:
km=1000;
z=linspace(-6000*km,0,600);
% the final time of the model:
year=3600*24*365;
time_end=20000e6*year; 
% time step:
dt=10^3*year;
% time axis for plotting:
time_axis=dt:dt:time_end;
% resolution in depth:
dz=(z(2)-z(1));
% number of time steps, and of grid points:
Nt=ceil(time_end/dt);
Nz=length(z);
% the depth where convection starts:
z_conv=-100*km; 
% diffusivities (m^2/s) are specified between temperature levels, such
% that D(k) is between levels k and k+1:
z_D(1:Nz-1)=(z(1:Nz-1)+z(2:Nz))/2;
% Kelvin diffusivity is constant in z:
D_K=1.e-6*ones(Nz-1,1);
% Perry diffusivity is larger in deeper earth by specified factor:
factor=700;
D_P=1.e-6*(factor*ones(Nz-1,1)-0.5*(factor-1)*(1+tanh((z_D-z_conv)/(3*dz))')); 
% plot diffusion coeffs for debugging purposes:
% figure(9);clf;plot(z_D,D_P,'r');hold on;plot(z_D,D_K,'b');legend('P','K');xlabel('z');return
% diffusion term= diffusivity times 2nd derivative WRT depth of
% temperature:
del2_temp=zeros(Nz,1);
% initial conditions for Perry solution; allow for three time levels
% for time stepping:
temp_P=4000*ones(Nz,3); 
% prescribed surface temperature:
temp_P(Nz,1:3)=273; 
% initial conditions for Kelvin's solution; allow for three time
% levels for time stepping:
temp_K=4000*ones(Nz,3); 
% prescribed surface temperature:
temp_K(Nz,1:3)=273; 

%% Integrate both cases in time:
n=1;
for t=0:dt:time_end
  n=n+1;
  %% (1) Perry's solution:
  %% ----------------------

  % center diffrencing of temperature:
  del2_temp(2:Nz-1)=( ...
    D_P(2:Nz-1).*(temp_P(3:Nz,1)  -temp_P(2:Nz-1,1)) - ...
    D_P(1:Nz-2).*(temp_P(2:Nz-1,1)-temp_P(1:Nz-2,1)) ...
    )/dz^2;
  % lower boundary (using no flux condition):
  del2_temp(1)=D_P(1)*(temp_P(2,1)-temp_P(1,1))/dz^2;
  % upper boundary (prescribed temperature):
  del2_temp(Nz)=0;
  
  % forwarding temp in time:
  temp_P(1:Nz-1,3)=temp_P(1:Nz-1,2)+del2_temp(1:Nz-1)*dt;

  % update time levels before next time step:
  temp_P(1:Nz-1,1)=temp_P(1:Nz-1,2);
  temp_P(1:Nz-1,2)=temp_P(1:Nz-1,3);
  
  % save T_z at z=0:
  grad_temp0_P(n)=(temp_P(Nz,3)-temp_P(Nz-1,3))/dz; 
  
  %% (2) Kelvin's solution:
  %% ----------------------
  
  % center diffrencing of temperature:
  del2_temp(2:Nz-1)=( ...
    D_K(2:Nz-1).*(temp_K(3:Nz,1)  -temp_K(2:Nz-1,1)) - ...
    D_K(1:Nz-2).*(temp_K(2:Nz-1,1)-temp_K(1:Nz-2,1)) ...
    )/dz^2;
  % lower boundary (using no flux condition):
  del2_temp(1)=D_K(1)*(temp_K(2,1)-temp_K(1,1))/dz^2;
  % upper boundary (prescribed temperature):
  del2_temp(Nz)=0;
  
  % forwarding temp in time:
  temp_K(1:Nz-1,3)=temp_K(1:Nz-1,2)+del2_temp(1:Nz-1)*dt;

  % update time levels before next time step:
  temp_K(1:Nz-1,1)=temp_K(1:Nz-1,2);
  temp_K(1:Nz-1,2)=temp_K(1:Nz-1,3);
  
  % save T_z at z=0:
  grad_temp0_K(n)=(temp_K(Nz,3)-temp_K(Nz-1,3))/dz; 

  
  %% Plot:
  %% -----
  %% set frequency of adding frames to animation:
  if t<4e9*year;
    nplot=20000;
  else
    nplot=200000;
  end
  %% add frame:
  if n==2 || floor(n/nplot)*nplot==n
    %% Kelvin:
    figure(1);clf
    pss=subplot(1,2,1);
    plot(temp_K(:,3),z/1e6,'b');
    set(gca,'fontsize',12);
    sdss=sprintf('Kelvin, T_z= %3.2g',abs(1.e3*grad_temp0_K(n)));
    title(sdss,'fontsize',16,'fontweight','bold');
    ylabel('Depth [10^3 km]','fontsize',16,'fontweight','bold');
    xlabel('Temperature [K]','fontsize',16,'fontweight','bold');
    xlim([0 4000]);
    
    %% time indicator on plot:
    text(3900,0.3,sprintf('t=%4.2fByr',t/(1.e9*year)) ...
         ,'fontsize',16,'fontweight','bold');
    
    %% Perry:
    pss=subplot(1,2,2);
    plot(temp_P(:,3),z/1e6,'r');
    set(gca,'fontsize',12);
    sdss=sprintf('Perry, T_z= %3.2g',abs(1.e3*grad_temp0_P(n)));
    title(sdss,'fontsize',16,'fontweight','bold');
    ylabel('Depth [10^3 km]','fontsize',16,'fontweight','bold');
    xlabel('Temperature [K]','fontsize',16,'fontweight','bold');
    xlim([0 4000]);

    pause(0.01);

    filename = ['Temperature_animation.gif'];
    %% write animation frames to file:
    if 1
      frame = getframe(1);
      im = frame2im(frame);
      [imind,cm] = rgb2ind(im,256);
      if n==2
        imwrite(imind,cm,filename,'gif', 'Loopcount',5,'DelayTime',0);
      else
        imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',0.1);
      end
    end % write animation to file
  
  end % plot
      
end % time loop

%return
%% Plot T_z at the surface as function of time for both solutions:
nfig=figure(2);clf; 
plot(time_axis(2:Nt)/(1e9*year),-grad_temp0_K(2:Nt)*1e3,'b','linewidth',2);
hold on;
plot(time_axis(2:Nt)/(1e9*year),-grad_temp0_P(2:Nt)*1e3,'r','linewidth',2);
ylim([0 200]);
set(gca,'fontsize',12);
ylabel('T_z [^\circ{}C/km]','fontsize',16,'fontweight','bold');
xlabel('Time [Byr]','fontsize',16,'fontweight','bold');
hl=legend('Kelvin - diffusion','Perry - diffusion and convection');
set(hl,'fontsize',16);
print (nfig,'-dpdf', 'T_z-timeseries-at-the-surface');
