%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program calculate the the contaminant transport equation with decay using 
% different schemes (explicit if e=0, implicit id e=1 and Crank-Nicholson if e=0.5).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
%
% define the different parameters
%
tmax=5;			    % simulation time - [day]
dx=0.5;				% spatial node - [m]
L=30;				% length of the domain - [m]
x=[0:dx:L]';        % length of the domain - [m]
V=2;				% velocity of the flow - [m/day]
C=0.1;				% Courant number
dt=C*dx/V;			% time space - [day]
t = 0:dt:tmax;      % time vector - [day]
P=10;				% Peclet number
D=V*dx/P;			% Diffusion coefficient - [m^2/day]
landa=0.5;			% Decay constant - [1/day]
Q =(D*landa)/(V^2); % Damkohler number
i=(L/dx)+1;			% size of the matrix
e=0.5;				% nature of the scheme (0:explicit,1:implicit,
				    % 0.5:Crank-Nicolson)
w=0.5;				% weighing parameter in space
c0=0;				% initial concentration at t=0 and x=0
c1=1;				% concentration of the contaminant at x=0 at all time
Flux = 0;           % Boundary flux at x = L
c=[c1;ones(i-1,1)*c0];	% initial concentration of the chemical reactant at t=0 for all x
ct=[c1;ones(i-1,1)*c0];	% initial concentration of the chemical tracer at t=0 for all x
%
% Determine the parameters and write them (diagonal terms) into the matrices M
% at time n+1 (Mn) and time n (Mo) => case for a chemical reactant
%
m=D/(dx^2);
n=-V/dx;

A=e*(m-n*w);
B=-1/dt-e*(2*m+n*(2*w-1)+landa);
E=e*(m+n*(1-w));
F=(e-1)*(m-n*w);
G=-1/dt+(1-e)*(2*m+n*(2*w-1)+landa);
H=(e-1)*(m+n*(1-w));

d1n=ones(i-1,1)*A;
d2n=ones(i,1)*B;
d3n=ones(i-1,1)*E;

Mn=diag(d2n,0)+diag(d1n,-1)+diag(d3n,+1);   % Building the matrix M at time n+1
Mn(1,:)=[1 zeros(1,i-1)];
Mn(i,:)=[zeros(1,i-2) -1 1];
 
d1o=ones(i-1,1)*F;
d2o=ones(i,1)*G;
d3o=ones(i-1,1)*H;  
    
Mo=diag(d2o,0)+diag(d1o,-1)+diag(d3o,+1);   % Building the matrix M a time n
Mo(1,:)=[1 zeros(1,i-1)];
Mo(i,:)=[zeros(1,i-2) -1 1];
%
% Determine the parameters and write them (diagonal terms) into the matrices M
% at time n+1 (Mn) and time n (Mo) => case for a chemical tracer
B1=-1/dt-e*(2*m+n*(2*w-1));
d2n=ones(i,1)*B1;
Mnt=diag(d2n,0)+diag(d1n,-1)+diag(d3n,+1);   % Building the matrix M at time n+1
Mnt(1,:)=[1 zeros(1,i-1)];
Mnt(i,:)=[zeros(1,i-2) -1 1];

G1=-1/dt+(1-e)*(2*m+n*(2*w-1));
d2o=ones(i,1)*G1;
Mot=diag(d2o,0)+diag(d1o,-1)+diag(d3o,+1);   % Building the matrix M a time n
Mot(1,:)=[1 zeros(1,i-1)];
Mot(i,:)=[zeros(1,i-2) -1 1];
%
% Calculate analytical solution for reactant and tracer
%
U = (V^2+4*D*landa)^0.5;
Ut = V;
%
% implement scheme for time and space and plot real time evolution of 
% the chemical reactant and the chemical tracer
%
for q=1:length(t)
    cn = inv(Mn)*(Mo*c);                      % concentration of the chemical reactant
    c = cn;
    
    ctn = inv(Mnt)*(Mot*ct);                  % concentration of the chemical reactant
    ct = ctn;
    
    A11 = 1/2*exp((V-U)*x/(2*D)).*erfc((x-U*t(q))/(2*(D*t(q))^0.5));
    A12 = 1/2*exp((V+U)*x/(2*D)).*erfc((x+U*t(q))/(2*(D*t(q))^0.5));
    Can = c1*(A11+A12);                   % analytical solution for the chemical reactant

    A21 = 1/2*exp((V-Ut)*x/(2*D)).*erfc((x-Ut*t(q))/(2*(D*t(q))^0.5));
    A22 = 1/2*exp((V+Ut)*x/(2*D)).*erfc((x+Ut*t(q))/(2*(D*t(q))^0.5));
    Ctan = c1*(A21+A22);
    
    plot(x,c,'r',x,Can,'b',x,ct,'r--',x,Ctan,'b--')
%    plot(x,c,'r',x,ct,'b--')
    title('Crank-Nicholson scheme for the transport of a chemical reactant and a tracer')
    xlabel('distance - [m]')
    ylabel('Conc. - [M]')
    text(17,0.8,'C = 0.1 ; P = 10.0 ; Q = 0.025')
    legend('reactant: numerical','reactant: analytical','tracer: numerical','tracer: analytical')
%    legend('reactant','tracer')
    pause(0.1)
end