% CH4_PROF an m-file to solve a boundary layer problem in anoxic
%	sediments. Problem taken from Martens and Berner (1977),
%	L&O, 22:10-25. The numerical integration is in addition
%	to the analytical solution in the paper.
%

% Set up some parameters

Ds=2e-6;	% Whole sediment diffusivity (cm^2 s^-1)
w=9.51e-9;	% Rate of burial (cm s^-1)
k1=7.61e-9; 	% First order rate constant (s^-1)
%k1=8e-9;
%w=0;
%k1=0;

% Set the grid spacing

m=35;
delx=175/m;

% Form the coefficients of the interior-point equations

A=(Ds/delx^2)*ones(1,m);
B=(w/delx+2*Ds/delx^2+k1)*ones(1,m);
C=(w/delx+Ds/delx^2)*ones(1,m);
D=0*ones(1,m);

% Set the boundary conditions

E1=0;
F1=0;
Wm=2600;

% Call the tridiagonal algorithm

[E,F,W]=tridiag(A,B,C,D,m,E1,F1,Wm);

% Generate a depth vector

z=0:delx:175-delx;

%Load the data

load ch4.dat;

% Plot the results

hold off;
plot(W,-z)
hold on;
plot(ch4(:,2),-ch4(:,1),'ro')
ylabel('Depth (cm)')
xlabel(['Methane (' setstr(181) 'M)'])