% heat1D.m
%
% Solves the 1D heat equation with an explicit finite difference scheme

clear

%Physical parameters
L       =   1e6;        %   length of modeled domain     [m]     
Tini    =   2300;       %   initial Earth temperature    [K]
Tsurf   =   300;        %   surface temperature          [K]
kappa   =   1e-6;       %   thermal diffusivity of rocks [m2/s]
day     =   3600*24;	%   # seconds per day
year    =   day*365;	%   # days per year
dt      =   1e5*year;	%   timestep                     [s]

% Numerical parameters
nz      =   101;	%   number of gridpoints in z-direction
nt      =   1000;	%   number of timesteps to compute
dz      =   L/(nz-1);	%   grid spacing "h"
z       =   0:dz:L;	%   grid

% Setup initial temperature profile
T       =   ones(size(z))*Tini;
T(1)    =   Tsurf;

% Setup time array to plot the thermal gradient
ti      =   zeros(1,nt);
Tgrad   =   zeros(1,nt);	% numerical solution
Tana    =   zeros(1,nt);	% analytical solution

time    =   0;
for n=1:nt		% Timestep loop

    % Compute new temperature 
    Tnew   =    zeros(1,nz);
    for i=2:nz-1
        Tnew(i) =  T(i) + kappa*dt*(T(i+1)+T(i-1)-2*T(i))/(dz*dz);
    end

    % Set boundary conditions
    Tnew(1)     =   T(1);
    Tnew(nz)    =   T(nz);
    
    % Update temperature and time
    T           =   Tnew;
    time        =   time+dt;
        
    % compute thermal gradient in K/km
    ti(n)       =   time;
    Tgrad(n)    =   1000*(T(2)-T(1))/dz;
    Tana(n)     =   1000*(Tini-Tsurf)/sqrt(pi*kappa*time);
    
    % print some progress information
    fprintf('%g out of %g: Tgrad = %f  Tana = %f  T(2) = %f\n', ...
             n,nt,Tgrad(n),Tana(n),T(2));
    
    % Plot solution
    figure(1), clf
    subplot(1,2,1)
    plot(T,-z/1000)
    grid
    xlabel('temperature [^oK]')
    ylabel('z [km]')
    title(['temperature T(z)'])
    subplot(1,2,2)
    plot(ti/(1e6*year),Tgrad)
    axis([0 100 0 100])
    xlabel('time [Myr]')
    ylabel('\DeltaT/\Deltaz_{num} [K km^{-1}]')
    title(['geothermal gradient'])
    
    drawnow
end
