%
% Example 4.2: Phosphorous Removal Process
%
% Bob Newell, March 1996
%
%-----------------------------------------------------------
%
function dxdt = CP_Removal( t, x )
global p
%
% x = [ 1=V 2=Sf 3=Sa 4=So 5=Spo 6=Xpha 7=Xpp 8=Xpao 9=Xh ];
% p = [ 1=muhath 2=Kf 3=Ka 4=Koh 5=Kp1 6=Yh 7=bh 8=fp
%       9=a 10=b 11=So_sat
%       12=muhatpao 13=Kpha 14=qfe 15=Kfe 16=qpha 17=Kpp
%       18=qpp 19=Kmax 20=Kipp 
%       21=bpha 22=bpp 23=bpao 24=Ypha 25=Ypo 26=ipbm
%       27=vp 28=airflow 29=Kp2 ];
%
%-----------------------------------------------------------
%
  muhetero = p(1) * x(4) * x(5) * x(9) ...
           / ( x(2) + x(3) ) ...
           / ( p(4) + x(4) ) / ( p(29) + x(5) );
  muheterof = muhetero * x(2) * x(2) / ( p(2) + x(2) );
  muheteroa = muhetero * x(3) * x(3) / ( p(3) + x(3) );
  mupao = p(12) * x(4) * x(5) * x(6) ...
           / ( p(4) + x(4) ) / ( p(29) + x(5) ) ...
           / ( p(13) + x(6) / x(8) );
  kla = p(9) * ( 1 - exp( - p(28) / p(10) ) );
  ferment = p(14) * p(4) * x(2) * x(9) ...
           / ( p(4) + x(4) ) / ( p(15) + x(2) );
  prelease = p(16) * x(3) * x(7) ...
           / ( p(3) + x(3) ) / ( p(17) + x(7) / x(8) );
  puptake = p(18) * x(4) * x(5) * x(6) ...
           * ( p(19) - x(7) / x(8) ) ...
           / ( p(4) + x(4) ) / ( p(5) + x(5) ) ...
           / ( p(13) + x(6) / x(8) ) ...
           / ( p(20) + p(19) - x(7) / x(8) );
%
% x = [ 1=V 2=Sf 3=Sa 4=So 5=Spo 6=Xpha 7=Xpp 8=Xpao 9=Xh ];
%
  dxdt = zeros( 9, 1 );
  dxdt(1) = 0;
  dxdt(2) = - ferment - muheterof + ( 1 - p(8) ) * p(7) * x(9);
  dxdt(3) =  ferment - muheteroa - prelease + p(21) * x(6);
  dxdt(4) = - p(24) * puptake + ( 1 - 1 / p(6) ) * mupao ...
            + ( p(6) - 1 ) * ( muheterof + muheteroa ) / p(6) ...
            + kla * ( p(11) - x(4) );
  dxdt(5) = p(25) * prelease - puptake ...
              - p(26) * mupao + p(2) * x(7) ...
              + p(27) * p(23) * x(8);
  dxdt(6) = prelease - p(24) * puptake ...
              - mupao / p(6) - p(21) * x(6);
  dxdt(7) = - p(25) * prelease + puptake ...
              - p(22) * x(7);
  dxdt(8) = mupao - p(23) * x(8);
  dxdt(9) = muheterof + muheteroa - p(7) * x(9);
%-----------------------------------------------------------
% the end!   
