%* 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;                       % number of 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;                     % 1/2 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)';  

%** numerics

%* set up A matrix (left-hand side of model equations)
A=zeros(N+1, N+1);          

% boundary values
A(1,1)= 1;
A(N+1,N+1)= 1;

% interior rows
for i=2:N
    A(i,i-1)= -K*dt/dz^2;
    A(i,i)= 1 + 2*K*dt/dz^2;
    A(i,i+1)= -K*dt/dz^2;
end

%** let's go!

for n=2:M                  % time steps, n=1 is initial condition
    
    % set up b vector (right-hand side of model equations) using T of time n-1 
    b=zeros(N+1,1);
    b(1)=Tsurf(n);      % boundary condition
    b(N+1)=Tbase;       % boundary condition
    b(2:N)=T(2:N,n-1);

    %solve
    T(:,n)=A\b;
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, implicit 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)')