%
% Example 14.5: Activated Sludge Aerobic Zone
%
% Batch reactor formulation, C only, air flow constant 
%
% Bob Newell, May 1996
%
%-----------------------------------------------------------
% parameters
muhath =6.0;
Kla = 75;
Ko2 =0.2;
Kf = 4.0;
Yh = 0.63;
%
global p
p = [ muhath Kla Ko2 Kf Yh ];
%-----------------------------------------------------------
% initial states
Sf = 60;             % 20 - 250
So2 = 0;
Xh = 300;
Xs = 150;
%
xo = [ Sf So2 Xh Xs ];
%-----------------------------------------------------------
% operational section
to = clock;
days = 8/24;
[ t, x ] = ode45( '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) );
%-----------------------------------------------------------
% plot the results
figure
hax = plotyy( 24*t, x(:,[1 3 4]), 24*t, [x(:,2) our/60] );
axesyy( hax, [ 0 24*days -500 500 0 30 ], ...
   'Time in hours', ...
   'Carbon in mg/l', ...
   'DO in mg/l & scaled OUR', 2, ...
   [ 0 250 500 ], [ 0 5 10 15 ] );
%
h = text( 2, 300, 'heterotrophs' );
set( h, 'fontsize', 16 );
h = text( 2, 170, 'insoluble COD' );
set( h, 'fontsize', 16 );
h = text( 2, -50, 'soluble COD' );
set( h, 'fontsize', 16 );
h = text( 5, -220, 'DO' );
set( h, 'fontsize', 16 );
h = text( 5, -420, 'OUR' );
set( h, 'fontsize', 16 );
%
%-----------------------------------------------------------
% the end!   
