%
% Example 3.6: Simple Biological Reactor
%
% - constant volume
% - well mixed
% - Monod kinetics
%
% Bob Newell, March 1996
%
%-----------------------------------------------------------
% parameters
muhat = 5;
Kn = 20;
Yb = 0.67;
Kfilter = 400;
%
global p
p = [ muhat Kn Yb Kfilter ];
%-----------------------------------------------------------
% initial states
V = 1.2;
Sn = 29.5;
Xb = 349;
%
xo = [ V Sn Xb Sn Xb Xb ];
%-----------------------------------------------------------
% integrate the state equation
to = clock;
[ t, x ] = ode45( 'Simple_Biological_Reactor', [0 1.5], xo );
etime( clock, to )
%-----------------------------------------------------------
% evaluate the outputs
y = Simple_Biological_Reactor_Outputs( t, x ); 
%-----------------------------------------------------------
% plot the results
figure;
[ hax ] = plotyy( 24*t, [ y(:,6) x(:,6) ], 24*t, y(:,5) );
axesyy( hax, [ 0 24 -300 400 0 100 ], ...
              'Time in hours', ...
              'Active biomass in mg/l', ...
              'Nutrient in mg/l', 0 );
%
figure;
n = length(t); m = 15;
plot( 24*t(1:m:n), y(1:m:n,6), '*', 24*t, x(:,6), '-' );

set( gca, 'fontsize', 16 );
xlabel( 'Time in hours' );
ylabel( 'Active biomass in mg/l' );
%
print -f2 -deps2 fig12-03.eps
%
%-----------------------------------------------------------
% integrate the state equation
p(4) = 4000;
to = clock;
[ t, x ] = ode45( 'Simple_Biological_Reactor', [0 1.5], xo );
etime( clock, to )
%-----------------------------------------------------------
% evaluate the outputs
y = Simple_Biological_Reactor_Outputs( t, x ); 
%-----------------------------------------------------------
% plot the results
figure;
[ hax ] = plotyy( 24*t, [ y(:,6) x(:,6) ], 24*t, y(:,5) );
axesyy( hax, [ 0 24 -300 400 0 100 ], ...
              'Time in hours', ...
              'Active biomass in mg/l', ...
              'Nutrient in mg/l', 0 );
%
figure;
n = length(t); m = 15;
plot( 24*t(1:m:n), y(1:m:n,6), '*', 24*t, x(:,6), '-' );

set( gca, 'fontsize', 16 );
xlabel( 'Time in hours' );
ylabel( 'Active biomass in mg/l' );
%
print -f4 -deps2 fig12-04.eps
%
%-----------------------------------------------------------
% the end!   
