%
% Exercise Program - Simulation of RAS Pipe
%
% Bob Newell, May 1996
%
%--------------------------------------------------
%
% standard conditions
%
nsect = 20;
volume = 6.28 / nsect;
n = 120;
disp( 'RAS Pipe Simulation' );
disp( ' ' );
vp1 = 60;
disp( 'Valve position   = 60 percent' );
ss1 = 8000;
disp( 'Suspended solids = 8000 mg/l' );
%
% request new values
%
vp2 = input( 'Enter new valve position ( 0 - 100 ) > ' );
if isempty( vp2 ), vp2 = vp1; end;
if vp2 < 0, vp2 = 0; end;
if vp2 > 100, vp2 = 100; end;
ss2 = input( 'Enter new concentration ( 5000 - 10000 )  > ' );
if isempty( ss2 ), ss2 = ss1; end;
if ss2 < 5000, ss2 = 5000; end;
if ss2 > 10000, ss2 = 10000; end;
%
% initialise states and calculate flowrate
% assuming equal percentage valve
%
for j = 1:nsect+1, x(j) = ss1; end;
flow = 0.52 * ( 8000 / ss1 ) * ...
                30 ^ ( vp1 / 100 - 1 );
%
% dynamic calculations
%
t0 = clock;
if ss2 ~= ss1,
  disp( 'Patience, watch for a plot ...' );
end;
kk = 1;
delt = 1;
for i = 1:n,
  if i == 9,
    x(1) = ss2;
    flow = 0.52 * ( 8000 / ss2 ) * ...
                30 ^ ( vp2 / 100 - 1 );
    kk = fix( flow * 8 + 1 );
    delt = 1 / kk;
  end;
  f(i) = flow;
  if ss2 ~= ss1,
    for k = 1:kk,
      for j = nsect+1:-1:2,
        x(j) = x(j) + delt * flow * ...
                    ( x(j-1) - x(j) ) / volume;
      end;
    end;
  end;
  ss(i) = x(nsect+1);
end;
secs = etime( clock, t0 )
%
% plot the results
%
figure;
plot( (1:n)*kk*delt, [ f' ss'/10000 ] )
axis( [ 0 n*kk*delt 0 1 ] );
set( gca, 'fontsize', 16 );
legend( 'Flow', 'SS', -1 );
xlabel( 'Time in seconds' );
%
%--------------------------------------------------
% the end
