%
%  Clarifier Model from Diehl and Jeppsson
%
%      [ xt, svf ] = Clarifier_Model( nn, ssi, ...
%                         fv, sv, cht, tht, delt, tmax, ...
%                         area, c, pflag, verbose )
%
%  where  xt      = suspended solids time & height profiles
%         svf     = suspended solids vector [ ssf sso ssu ]
%         nn      = number of nodes
%         ssi     = initial suspended solids (length nn)
%         fv      = flowrate vector [ feed overflow underflow ]
%         sv      = suspended solids vector [ ssf 0 0 ]
%         cht     = clarifier height
%         tht     = thickener height
%         delt    = time increment
%         tmax    = end of calculation time
%         area    = cross-sectional area
%         c       = coefficient vector for settling velocity model
%         pflag   = plot time (times delt)
%         verbose = flux and solids printout (times delt)
%
% Bob Newell, April 1996
%
%-------------------------------------------------------------------
function [ xt, svf ] = Clarifier_Model( nn, ssi, ...
                         fv, sv, cht, tht, delt, tmax, ...
                         area, c, pflag, verbose )
%
global coef
global bulkslope
%
% set up solution vectors
%
t = 0;
x = ssi(:);
xt = [ t x' ];
delh = ( cht + tht ) / ( nn - 2 );
nf = round( cht/delh + 1.5 );
f = zeros( nn-nf+1, 1 );
g = zeros( nf, 1 );
FG = zeros( nn, 1 );
%
% calculate xm and xM
%
disp( 'Searching for limiting flux region...' );
coef = c;
bulkslope = -1;
% find peak of flux curve without bulk flux
xm = fmin( 'Flux_Model', 0, 20000 );
% find minimum in flux curve with bulk flow
bulkslope = fv(3) / area;
xM = fmin( 'Flux_Model', xm, 30000 );
if xM == 30000,
  disp( 'Error - limiting flux region not found' );
  return;
end;
% now find bottom end of limiting flux region
fM = Flux_Model( xM );
while Flux_Model( xm ) > fM,
  xm = xm - 1;
end;
bulkslope = 0;
disp( [ '  from ' num2str( fix( xm ) ) ...
                       ' to ' num2str( fix( xM ) ) ] );
%
% plot initial solids profile
%
if pflag > 0,
  figure
  plot( x, nn-1:-1:0, '--' );
  xlabel( 'Solids' );
  ylabel( 'Height' );
  hold on;
end;
%
% integration of time
%
n = 0;
while t < tmax,
  t = t + delt;
  n = n + 1;
  [ fv, sv ] = Clarifier_Disturbance( t, fv, sv );
%
% calculate fluxes calculations
%
  vsx = Flux_Model( x );
%
% clarification zone flux calculations
%
  g(1:nf) = vsx(1:nf) - fv(2) * x(1:nf) / area;
%
% thickening zone flux calculations
%
  f(1:nn-nf+1) = vsx(nf:nn) + fv(3) * x(nf:nn) / area;
%
% total flux calculations in clarification zone
%
  for i = 2:nf,
    if x(i-1) <= x(i),
      FG(i-1) = min( g(i-1), g(i) );
    else,
      FG(i-1) = max( g(i), g(i-1) );
    end;
  end;
%
% total flux calculations in thickening zone
%
  for i = nf+1:nn,
    if x(i-1) <= x(i),
      FG(i-1) = min( f(i-nf), f(i-nf+1) );
    else,
      FG(i-1) = max( f(i-nf+1), f(i-nf) );
    end;
  end;
%
% time update of solids profile
%
  for i = 2:nn-1,
    x(i) = x(i) + delt * ( FG(i-1) - FG(i) ) / delh; 
  end;
  x(nf) = x(nf) + delt * fv(1) * sv(1) / area / delh;
%
% boundary values
%
  if g(2) <= 0,
    x(1) = x(2);
  else,
    x(1) = 0;
  end;
  if x(nn-1) > xm & x(nn-1) < xM,
    x(nn) = xM;
  else,
    x(nn) = x(nn-1);
  end;
%
% update plot of solids profile
%
  xt = [ xt; [ t x' ] ];
  fprintf( 1, '\nTime %f/%f', t, tmax );
  if verbose > 0,
    if verbose*(fix(n/verbose)+1) == n+1,
      fprintf( 1, '\nFlux (g) profile is:' );
      fprintf( 1, '\n%f  %f  %f  %f  %f', g );
      fprintf( 1, '\nFlux (f) profile is:' );
      fprintf( 1, '\n%f  %f  %f  %f  %f', f );
      fprintf( 1, '\nFlux (FG) profile is:' );
      fprintf( 1, '\n%f  %f  %f  %f  %f', FG );
      fprintf( 1, '\nSolids profile is:' );
      fprintf( 1, '\n%f  %f  %f  %f  %f', x );
      fprintf( 1, '\n' );
      pause;
    end;
  end;
  if pflag > 0,
    if pflag*(fix(n/pflag)+1) == n+1,
      plot( x, nn-1:-1:0, '--' );
      pause
    end;
  end;
%
% calculate effluent suspended solids
%
  sv(2) = x(1) - vsx(1) * area / fv(2);
  sv(3) = x(nn) + vsx(nn) * area / fv(3);
%
end;
%
% update plot of solids profile
%
% disp( FG );
if pflag > 0,
  plot( x, nn-1:-1:0, 'm-' );
  hold off;
end;
%
fprintf( 1, '\nTime is %f, Final solids profile is:\n', t );
fprintf( 1, '%10d%10d%10d%10d%10d\n', fix( x ) );
fprintf( 1, 'Overflow is  %d\n', fix( sv(2) ) );
fprintf( 1, 'Underflow is %d\n', fix( sv(3) ) );
%
xt = [ xt; [ t x' ] ];
svf = sv;
return;
%----------------------------------------------------------------
