% ps08_2.m is  an m-file to solve problem #2 from problem set 
% #08 which entails modeling Radon 222 in the sediments.
%
% Modif'd 2002:11:20 D. Glover for 2002

% Set up some parameters (given in the problem)

Ds=2.66e-6;	% Whole sediment diffusivity (cm^2 s^-1)
w=9.6e-10;	% Rate of burial (cm s^-1)
k1=2.098e-6; 	% Decay constant of 222Rn (s^-1)
k2=1.373e-11;   % Decay constant of 226Ra (s^-1)
Ra=9.614e7;	% 226Ra concentration (constant over sed column)
Rn=283.1;
%w=0;
%k1=0;

% Set the grid spacing

m=90;
delx=30/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=(Ra*k2)*ones(1,m);
Sec=Ra*k2*ones(1,m);

% Set the boundary conditions

E1=0;
F1=Rn;
Wm=Ra*k2/k1;

% Call the tridiagonal algorithm

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

% Generate a depth vector

z=0:delx:30-delx;

% Plot the results
figure(1);clf;
hold off;
plot(W*k1,z)
hold on;
plot(Sec,z,'--')
axis('ij')
ylabel('Depth (cm)')
xlabel(['222 Radon Activity  (dps/cm^3)'])
title('^{222}Rn Activity vs. Depth (No Bio-irrigation)')
Flux=trapz(z,(W*k1-Wm*k1))
text(5.5e-4,10,(sprintf('^{222}Rn Flux = %d (atoms cm^{-2} s^{-1})',Flux)))



% Now for part b (bio-irrigation)
% Set up some parameters (given in the problem)

Ds=2.66e-6;   % Whole sediment diffusivity (cm^2 s^-1)
w=9.6e-10;    % Rate of burial (cm s^-1)
VB=9.3e-6;    % Bio-irrigation (actually 6.51e-6 (cm s^1) accounting for n=.7. 
k1=2.098e-6;  % Decay constant of 222Rn (s^-1)
k2=1.373e-11; % Decay constant of 226Ra (s^-1)
Ra=9.614e7;   % 226Ra concentration (constant over sed column)
Rn=283.1;
%w=0;
%k1=0;

% Set the grid spacing

m=90;
delx=30/m;

% Form the coeeficients of the interior-point equations

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

% Set the boundary conditions

E1=0;
F1=Rn;
Wm=Ra*k2/k1;

% Call the tridiagonal algorithm

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

% Generate a depth vector

z=0:delx:30-delx;

% Plot the results
figure(2);clf;
hold off;
plot(W*k1,z)
hold on;
plot(Sec,z,'--')
axis('ij')
ylabel('Depth (cm)')
xlabel(['222 Radon Activity  (dps/cm^3)'])
title('^{222}Rn Activity vs. Depth (With Bio-irrigation)')
Flux=trapz(z,(W*k1-Wm*k1))
text(6e-4,15,(sprintf('^{222}Rn Flux = %d (atoms cm^{-2} s^{-1})',Flux)))
VBI=VB*.7*8.64e4;
text(6e-4,20,(sprintf('V_{BI} = %d (cm d^{-1})',VBI)))

disp('The deficit (i.e. flux) in part (a) was 8.3e-4 atoms cm-2 s-1,');
disp('considerably smaller than the sediment core measurements estimate.');
disp('Yes, with a relatively small amount of bioirrigation, the flux ');
disp('of 222Rn can easily be raised to the observed 3.5e-3 atoms cm-2 s-1.');
disp('The bioirrigation necessary was only 0.56 cm d-1, compared to ');
disp('measurements on the Washington continental shelf this is only a');
disp('quarter of what is observed. Since the Bering Sea shelf is farther');
disp('north, this lower bioirrigation is acceptable. VBI is divided by');
disp('porosity because only the water is irrigated through the sediment, ');
disp('the sediment itself stays behind.');