%
% Example 4.1: Nitrogen Removal Process
%
% Bob Newell, March 1996
%
%-----------------------------------------------------------
%
function dxdt = CN_Removal( t, x )
global p
%
% x = [ 1=V 2=Ss 3=So 4=Snh 5=Sno 6=Xh 7=Xa ];
% y = [ 1=v_a ];
% p = [ 1=muhath 2=Ks 3=Koh 4=Yh 5=bh 6=fp
%       7=a 8=b 9=So_sat
%       10=muhata 11=Knh 12=Koa 13=Ya 14=ba 
%       15=Kno 16=etag 17=ixb 18=ixp 19=airflow ];
%
%-----------------------------------------------------------
%
  muhetero = p(1) * x(2) * x(3) * x(6) / ( p(2) + x(2) ) ...
                         / ( p(3) + x(3) );
  muanoxic = p(1) * x(2) * p(3) * x(5) * p(16) * x(6) ...
                / ( p(2) + x(2) ) / ( p(15) + x(5) ) ...
                         / ( p(3) + x(3) );
  muauto = p(10) * x(4) * x(3) * x(7) ...
                / ( p(11) + x(4) ) / ( p(12) + x(3) );
  kla = p(7) * ( 1 - exp( - p(19) / p(8) ) );
%
% x = [ 1=V 2=Ss 3=So 4=Snh 5=Sno 6=Xh 7=Xa ];
%
  dxdt = zeros( 7, 1 );
  dxdt(1) = 0;
  dxdt(2) = - muhetero / p(4) - muanoxic / p(4) ...
                + p(5) * x(6) * ( 1 - p(6) ) ...
                + p(14) * x(7) * ( 1 - p(6) );
  dxdt(3) = muhetero * ( p(4) - 1 ) / p(4) ...
                + muauto * ( p(13) - 4.57 ) / p(13) ...
                + kla * ( p(9) - x(3) );
  dxdt(4) = - muhetero * p(17) - muanoxic * p(17) ...
                - muauto * ( p(17) + 1 / p(13) ) ...
                + ( p(17) - p(6) * p(18) ) ...
                * ( p(5) * x(6) + p(14) * x(7) );
  dxdt(5) = muanoxic * ( p(4) - 1 ) / 2.86 / p(4) ...
                + muauto / p(13);
  dxdt(6) = muhetero + muanoxic - p(5) * x(6);
  dxdt(7) = muauto - p(14) * x(7);
%-----------------------------------------------------------
% the end!   
