


% This is a .m matlab file run it within matlab





disp(' -------------------------------------------------- ')
disp('')
disp('')
disp('                 Diffusion-Advection-Reaction      ')
disp('')
disp('')
disp('  ')
disp('')


pause off


% initialize constants

k = 3.e4;     	  % Eddie mixing in km^2/a
v = 1500.;        % velocity km/a
p0 = 0.01;           % production rate Mol/a
r = 2e-3;        % removal rate constant I guess Mol/a

% space-time discretization (and misc)


nt = 30000;             % number of time steps
nx = 100;              % space dimension

phi=2*pi/(nx-1)*(0:nx-1);
xc=cos(phi);
yc=sin(phi);


dt = 0.01;   	   	% years
dx = 1000;                % km 
isnap = 50;             % snapshot sampling rate


% initialize space variables
c  = zeros([1 nx]);              % set c to zero
dc = zeros([1 nx]);              % dc/dx
p=(1:nx);p=p*0;p(nx/2)=p0;


% Time extrapolation and plotting

disp(' Beginning time extrapolation .....')

for it = 1:nt,

% calculate the rhs

%dc=k*(csh(c,1)-2*c+csh(c,-1))/dx^2+v*(csh(c,1)-csh(c,-1))/(2*dx)+p-r*c; % centered (velocity term)
%dc=k*(csh(c,1)-2*c+csh(c,-1))/dx^2+v*(csh(c,0)-csh(c,-1))/(dx)+p-r*c; % backward(velocity term)
dc=k*(csh(c,1)-2*c+csh(c,-1))/dx^2+v*(csh(c,1)-csh(c,0))/(dx)+p-r*c; % forward (velocity term)  - THE WINNER

% extrapolate to next time step

c=c + dc * dt;

% spatial boundary condition
c(1)=c(nx);

% visualization

if rem(it,isnap) == 0,
%plot((1:nx)*dx,c,'b-'),
plot3(xc,yc,c,'b-',xc,yc,c*0.,'r-')
title(sprintf(' Concentration in year %g', it*dt))
grid
drawnow

end
end 