%
% Example 4.1: Nitrogen Removal Process
%
% - constant volume
% - well mixed
% - Monod kinetics
% - batch configuration
%
% Bob Newell, March 1996
%
%-----------------------------------------------------------
% parameters
muhath = 5;
Ks = 20;
Koh = 0.2;
Yh = 0.63;
bh = 0.4;
fp = 0.1;
a = 166;
b = 16;
So_sat = 10;
muhata = 1;
Knh = 1;
Koa = 0.5;
Ya = 0.24;
ba = 0.15;
Kno = 0.5;
etag = 0.8;
ixb = 0.1;
ixp = 0.1;
air = 50;
%
global p
p = [ muhath Ks Koh Yh bh fp a b So_sat ...
      muhata Knh Koa Ya ba Kno etag ixb ixp air ];
%-----------------------------------------------------------
% initial states
V = 12.5;
Ss = 150;
So = 3;
Snh = 40;
Sno = 0;
Xh = 400;
Xa = 40;
%
xo = [ V Ss So Snh Sno Xh Xa ];
%-----------------------------------------------------------
% integrate the state equation
to = clock;
days1 = 0.333333;
[ t, x ] = ode45( 'CN_Removal', [ 0 days1 ], xo );
secs = etime( clock, to )
%-----------------------------------------------------------
% plot the results
figure;
hax = plotyy( 24*t, x(:,[6 7]), 24*t, [x(:,2)/10 x(:,[3 4 5])] );
axesyy( hax, [ 0 24*days1+.01 -800 800 0 80 ], ...
   'Time in hours', ...
   'Active biomass in mg/l', ...
   'Nutrients in mg/l', 2, ...
   [ 0 200 400 600 ], [ 0 20 40 ] );
%
stext( 2, 550, '{\times \16 \it X_H}' );
stext( 2, 100, '{\times \16 \it X_A}' );
%
stext( 0.5, -550, '{\times \16 \it S_S /10}' );
stext( 6.5, -540, '{\times \16 \it S_O}' );
stext( 1.8, -280, '{\times \16 \it S_{NH4}}' );
stext( 6.5, -140, '{\times \16 \it S_{NOX}}' );
%
printsto -f1 -deps2 ex4-1a.eps
%
%-----------------------------------------------------------
xo = x(length(x),:);
xo(2) = 400;
xo(4) = 35;
p(19) = 0;
%-----------------------------------------------------------
% integrate the state equation
to = clock;
days2 = 0.5;
[ t, x ] = ode45( 'CN_Removal', [ days1 days2 ], xo );
secs = etime( clock, to )
%-----------------------------------------------------------
% plot the results
figure;
hax = plotyy( 24*t, x(:,[6 7]), 24*t, [x(:,2)/10 x(:,[3 5])] );
axesyy( hax, [ 24*days1 24*days2 -800 800 0 80 ], ...
   'Time in hours', ...
   'Active biomass in mg/l', ...
   'Nutrients in mg/l', 2, ...
   [ 0 200 400 600 ], [ 0 20 40 ] );
%
stext( 9, 620, '{\times \16 \it X_H}' );
stext( 9, 130, '{\times \16 \it X_A}' );
%
stext( 9.5, -300, '{\times \16 \it S_S /10}' );
stext( 8.1, -680, '{\times \16 \it S_O}' );
stext( 8.8, -600, '{\times \16 \it S_{NOX}}' );
%
printsto -f2 -deps2 ex4-1b.eps
%-----------------------------------------------------------
% the end!   
