%
% Example 4.2: Phosphorous Removal Process
%
% - constant volume
% - well mixed
% - Monod kinetics
% - batch configuration
%
% Bob Newell, March 1996
%
%-----------------------------------------------------------
% parameters
muhath = 5;
Kf = 4;
Ka = 4;
Koh = 0.2;
Kp1 = 0.2;
Yh = 0.63;
bh = 0.4;
fp = 0.1;
a = 166;
b = 16;
So_sat = 10;
muhatpao = 1;
Kpha = 0.01;
qfe = 3;
Kfe = 20;
qpha = 3;
Kpp = 0.01;
qpp = 1.5;
Kmax = 0.34;
Kipp = 0.02;
bpha = 0.2;
bpp = 0.2;
bpao = 0.2;
Ypha = 0.2;
Ypo = 0.4;
ipbm = 0.02;
vp = 0.01;
air = 0;
Kp2 = 0.01;
%
global p
p = [ muhath Kf Ka Koh Kp1 Yh bh fp a b So_sat ...
      muhatpao Kpha qfe Kfe qpha Kpp qpp Kmax Kipp ...
      bpha bpp bpao Ypha Ypo ipbm vp air Kp2 ];
%-----------------------------------------------------------
% initial states
V = 12.5;
Sf = 350;
Sa = 50;
So = 3;
Spo = 5;
Xpha = 0;
Xpp = 35;
Xpao = 700;
Xh = 150;
%
xo = [ V Sf Sa So Spo Xpha Xpp Xpao Xh ];
%-----------------------------------------------------------
% integrate the state equation
to = clock;
days1 = 1/24;
[ t, x ] = ode45( 'CP_Removal', [ 0 days1 ], xo );
secs = etime( clock, to )
%-----------------------------------------------------------
% plot the results
font = 20;
figure;
hax = plotyy( 24*t, [x(:,[9 8]) 10*x(:,6)]/1000, ...
   24*t, [x(:,2)/4 x(:,[3 4 5 7])] );
axesyy( hax, [ 0 24*days1 -1.000 1.000 0 200 ], ...
   'Time in hours', ...
   'Active biomass in g/l', ...
   'Nutrients in mg/l', 2, ...
   [ 0 0.5 1.0 ], [ 0 50 100 ] );
%
stext( 0.8, .250, '{\times \20 \it X_H}' );
stext( 0.2, .800, '{\times \20 \it X_{PAO}}' );
stext( 0.4, .450, '{\times \20 \it X_{PHA}*10}' );
%
stext( 0.4, -.250, '{\times \20 \it S_F /4}' );
stext( 0.1, -.500, '{\times \20 \it S_A}' );
stext( 0.8, -.630, '{\times \20 \it S_{PO4}}' );
stext( 0.05, -.750, '{\times \20 \it X_{PP}}' );
%
printsto -f1 -deps2 ex4-2a.eps
%
%-----------------------------------------------------------
xo = x(length(x),:);
%xo(2) = 260;
p(28) = 50;
%-----------------------------------------------------------
% integrate the state equation
to = clock;
days2 = 4/24;
[ t, x ] = ode45( 'CP_Removal', [ days1 days2 ], xo );
secs = etime( clock, to )
%-----------------------------------------------------------
% plot the results
figure;
hax = plotyy( 24*t, [x(:,[9 8]) 10*x(:,6)]/1000, ...
   24*t, [x(:,2)/4 x(:,[3 4 5 7])] );
axesyy( hax, [ 24*days1 24*days2 -1.000 1.000 0 200 ], ...
   'Time in hours', ...
   'Active biomass in g/l', ...
   'Nutrients in mg/l', 2, ...
   [ 0 0.5 1.0 ], [ 0 50 100 ] );
%
stext( 3, .310, '{\times \20 \it X_H}' );
stext( 2.4, .820, '{\times \20 \it X_{PAO}}' );
stext( 1.7, .450, '{\times \20 \it X_{PHA}*10}' );
%
stext( 2.8, -.200, '{\times \20 \it S_F /4}' );
stext( 2.5, -.870, '{\times \20 \it S_O}' );
stext( 1.2, -.890, '{\times \20 \it S_A}' );
stext( 1.8, -.800, '{\times \20 \it S_{PO4}}' );
stext( 3.5, -.835, '{\times \20 \it S_{PO4}}' );
stext( 1.6, -.500, '{\times \20 \it X_{PP}}' );
%
printsto -f2 -deps2 ex4-2b.eps
%
%-----------------------------------------------------------
% the end!   
