%
% Example 18.4: GMC of DO on Carbon Removal Process
%
% Bob Newell, April 1996
%
%-----------------------------------------------------------
%
function y = Carbon_Removal_Outputs( t, x )
global p
%
% x = [ 1=V 2=Sn 3=So 4=Xb 5=z 6=OURmeas ];
% y = [ 1=v_in 2=Ss_in 3=So_in 4=Xh_in 
%       5=v_out 6=Ss_out 7=So_out 8=Xh_out
%       9=v_a 10=v_r 11=v_w 12=OUR ];
% p = [ 1=muhat 2=Ks 3=Koh 4=Yh 5=bh 6=fp
%       7=a 8=b 9=So_sat
%       10=So_desired 11=K1 12=K2 13=Kp 14=Ki 15=law ];
%
%-----------------------------------------------------------
  [ n, m ] = size( x );
  if m == 1,
    x = x';
    n = 1;
  end;
  y = zeros( n, 2 );
%
  for i = 1:n,
    % influent conditions
    y(i,3) = 1;
    y(i,4) = 5;
    if t(i) < 1/24,
      y(i,1) = 50;
      y(i,2) = 150;
      y(i,9) = 6.244;
    elseif t(i) < 3.5/24,
      y(i,1) = 75;
      y(i,2) = 200;
      y(i,9) = 6.244;
    elseif t(i) < 6/24,
      y(i,1) = 25;
      y(i,2) = 100;
      y(i,9) = 6.244;
    else,
      y(i,1) = 50;
      y(i,2) = 150;
      y(i,9) = 6.244;
    end; 
    % sludge flows
    y(i,10) = 50;
    y(i,11) = 1.615;
    % constant volume
    y(i,5) = y(i,1);
    % well mixed
    y(i,6) = x(i,2);
    y(i,7) = x(i,3);
    y(i,8) = x(i,4);
    mu = p(1) * y(i,6) * y(i,7) * y(i,8) ...
	      / ( p(2) + y(i,6) ) / ( p(3) + y(i,7) );
    y(i,12) = - ( p(4) - 1 ) * mu / p(4);
    if p(15) < 2,
% this is a PI controller
      y(i,9) = y(i,9) + p(13) * ( p(10) - y(i,7) ) ...
	                + p(14) * x(i,5);
    elseif p(15) < 3,
% this is the GMC controller with OUR
      y(i,14) = p(11) * ( p(10) - y(i,7) );
      y(i,15) = p(12) * x(i,5);
      term1 = p(7) * ( p(9) - y(i,7) );
      term2 = ( y(i,1) * y(i,3) ...
        - ( y(i,1) + y(i,10) ) * y(i,7) ) / x(i,1) ...
    	+ term1 - x(i,6) - y(i,14) - y(i,15);
      y(i,9) = p(8) * log( term1 / term2 );
    else,
% no control at all
    end;
  end;
%-----------------------------------------------------------
% the end!   
