%* 1-D thermal diffusion model
clear

%** define constant(s)
K=2.8e-7;                  % thermal diffusivity of dry sand, m^2/s 
day=24*60*60;           % seconds in one day
year=365;               % days in one year 


%** model domain
zsurf=0;                        % surface coordinate, m
zbase=10;                       % base coordinate, m  "deep" 
dz=0.1;                         % vertical interval size, m
z=[zsurf:dz:zbase]';            % vertical coordinates, m
N=length(z)-1;                  % intervals in model domain

ti=0;               % initial time
tf=1;               % final time in years 
nD=0.1;             % time step size in days

dt=nD*day;          % time step size in seconds
M=tf*year/nD;          % number of time steps, days 


%* boundary conditions
%  start model at coldest part of temperature cycle, -1/2 pi
Tsurfmean=280;                % mean annual temperature, K
Tdiff=30;                     % annual temperature range, K
Tsurf=Tsurfmean + Tdiff/2*sin(linspace(-1/2*pi, (tf*2-1/2)*pi, M)); 
Tbase=Tsurfmean;  

%* variables
T=zeros(N+1,M);           %T history 

%* insert boundary and initial conditions into matrix T 
%  initial condition: use first Tsurf & monotonic trend to Tbase 
T(1,1:M)=Tsurf;
T(N+1,1:M)=Tbase;

T(:,1)=linspace(Tsurf(1),Tbase, N+1)';  

%** let's go!

for n=1:M-1                  % time steps, n=1 is initial condition
    for j=2:N                  % space steps, N+1 nodes
        T(j, n+1) = T(j,n) + K*dt/dz^2 * (T(j+1,n) - 2*T(j,n) + T(j-1,n));
    end
end



%** plot result

%* a family of curves representing T(z) at selected time steps
Tmin=min(min(T));
Tmax=max(max(T));

figure(3)
clf
axis([Tmin Tmax -zbase zsurf])      % use -z to plot in a perspective that makes physical sense
hold on
plot(T(:,2:30/nD:M-1), -z, 'c-')   % plot every 30th day
plot(T(:,1), -z, 'b-')              % plot initial condition
ylabel('depth (m)')
xlabel('T (K)')
title('temperature cycle in a soil layer, explicit numerical scheme')



%* a plot of the complete T(z,t) space
%figure(4)
%surf(z,[ti:M-1]*nD,T')
%shading interp
%xlabel('depth (m)')
%ylabel('day of year')
%zlabel('T (K)')