%
% Exercise 7: Carbon Removal Process
%
% Bob Newell, May 1996
%
%-----------------------------------------------------------
%
function dxdt = Carbon_Plant( t, x )
global p
%
% x = [ 1=Ss1 2=So1 3=Ss2 4=So2 5=Ss3 6=So3 ];
% y = [ 1=v_in 2=Ss_in 3=So_in 4=Xh_in
%       5=Ss1 6=Ss2 7=Ss3 8=So1 9=So2 10=So3 ];
% p = [ 1=muhat 2=Ks 3=Koh 4=Yh 5=bh 6=fp
%       7=a 8=b 9=So_sat
%       10=vol1 11=vol2 12=vol3
%       13=Kla1 14=Kla2 15=Kla3
%       16=flow 17=cod ];
%
%-----------------------------------------------------------
y = Carbon_Plant_Outputs( t, x );
  dxdt = zeros( 6, 1 );
%
  mu = p(1) * x(1) * x(2) * y(4) / ( p(2) + x(1) ) ...
                         / ( p(3) + x(2) );
  dxdt(1) = y(1) * ( y(2) - x(1) ) / p(10) ...
                - mu / p(4) + p(5) * y(4) * ( 1 - p(6) );
  dxdt(2) = y(1) * ( y(3) - x(2) ) / p(10) ...
                + mu * ( p(4) - 1 ) / p(4) ...
                + p(13) * ( p(9) - x(2) );
%
  mu = p(1) * x(3) * x(4) * y(4) / ( p(2) + x(3) ) ...
                         / ( p(3) + x(4) );
  dxdt(3) = y(1) * ( x(1) - x(3) ) / p(11) ...
                - mu / p(4) + p(5) * y(4) * ( 1 - p(6) );
  dxdt(4) = y(1) * ( x(2) - x(4) ) / p(11) ...
                + mu * ( p(4) - 1 ) / p(4) ...
                + p(14) * ( p(9) - x(4) );
%
  mu = p(1) * x(5) * x(6) * y(4) / ( p(2) + x(5) ) ...
                         / ( p(3) + x(6) );
  dxdt(5) = y(1) * ( x(3) - x(5) ) / p(12) ...
                - mu / p(4) + p(5) * y(4) * ( 1 - p(6) );
  dxdt(6) = y(1) * ( x(4) - x(6) ) / p(12) ...
                + mu * ( p(4) - 1 ) / p(4) ...
                + p(15) * ( p(9) - x(6) );
%
%-----------------------------------------------------------
% the end!   
