%
% Example 14.6: Activated Sludge Aerobic Zone
%
% Batch reactor formulation, C & N only, air flow constant 
%
% Bob Newell, May 1996
%
%-----------------------------------------------------------
% parameters
muhath =6.0;
Kla = 100;
Ko2 =0.2;
Kf = 4.0;
Yh = 0.63;
muhata = 1.0;
Ko22 = 0.5;
Knh4 = 1;
Ya = 0.24;
%
global p
p = [ muhath Kla Ko2 Kf Yh muhata Ko22 Knh4 Ya ];
%-----------------------------------------------------------
% initial states
Sf = 40;             % 20 - 250
So2 = 0;
Xh = 300;
Xs = 50;
Snh4 = 15;
Sno3 = 0;
Xa = 30;
%
xo = [ Sf So2 Xh Xs Snh4 Sno3 Xa ];
%-----------------------------------------------------------
% operational section
to = clock;
days = 8/24;
[ t, x ] = ode45( 'Nitrifying_Aerobic_Zone', [ 0 days ], xo );
secs = etime( clock,to)
%-----------------------------------------------------------
% calculate the OUR
our = - p(1) * x(:,1) .* x(:,2) .* x(:,3) ...
              * ( p(5) - 1 ) / p(5) ...
              ./ ( p(3) + x(:,2) ) ./ ( p(4) + x(:,1) );
nour = - p(6) * x(:,2) .* x(:,5) .* x(:,7) ...
              * ( p(9) - 4.57 ) / p(9) ...
              ./ ( p(7) + x(:,2) ) ./ ( p(8) + x(:,5) ) ;
%-----------------------------------------------------------
% plot the results
figure
hax = plotyy( 24*t, x(:,[1 4 5]), ...
              24*t, [x(:,2) our/100 nour/100] );
axesyy( hax, [ 0 24*days -60 60 0 20 ], ...
   'Time in hours', ...
   'Carbon and nitrogen in mg/l', ...
   'DO in mg/l and scaled OURs', 2, ...
   [ 0 20 40 60 ], [ 0 4 8 ] );
%
h = text( 0.8, -25, 'COUR' );
set( h, 'fontsize', 16 );
h = text( 3, -32, 'NOUR' );
set( h, 'fontsize', 16 );
h = text( 6, -28, 'DO' );
set( h, 'fontsize', 16 );
h = text( 1.5, -6, 'soluble COD' );
set( h, 'fontsize', 16 );
h = text( 4.5, 13, 'insoluble COD' );
set( h, 'fontsize', 16 );
h = text( 2, 40, 'ammonia' );
set( h, 'fontsize', 16 );
arrow( [2 36], [1.5 12] );
%
printsto -f1 -deps2 fig14-13.eps
%
%-----------------------------------------------------------
% the end!   
