%
% RTD Simulator - Tracer injection simulation
%
% Bob Newell, January 1996
%
function results = rtdcalc( objects )
[nobjects, n] = size( objects );
nsect = 50;
states = zeros( nobjects, nsect+2 );
%
% this is our integrator (if you can call it that?)
%
disp( 'Patience, computer at work ...' );
maxtime = 3;
delt = 0.01;
results = zeros( maxtime/delt+1, nobjects+1 );
n = 2;
for time = delt:delt:maxtime,
%
% calculate inlet flows (total and tracer)
%
  
  for i = 1:nobjects,
    states(i,1) = 0;
    states(i,2) = 0;
  end;
  for i = 1:nobjects,
    for j = 2:5,
      k = objects(i,2*j);
      if k > 0,
        states(k,1) = states(k,1) + objects(i,2*j+1);
        states(k,2) = states(k,2) + ...
                objects(i,2*j+1) * states(i,nsect+2);
      end;
    end;
  end;
%
% calculate next value of state(s)
%
  for i = 1:nobjects,
    if objects(i,1) == 1,
%
% tracer injection
%
      if objects(i,2) == 0,
        if time < delt * 2,
          states(i,nsect+2) = objects(i,3);
        else
          states(i,nsect+2) = 0;
        end;
      elseif objects(i,2) == 1,
        if time < delt * 50,
          states(i,nsect+2) = objects(i,3);
        else
          states(i,nsect+2) = 0;
        end;
      elseif objects(i,2) == 2,
        states(i,nsect+2) = objects(i,3);
      end;
%
    elseif objects(i,1) == 2,
%
% well-mixed tank
%
      dcdt(1) = ( states(i,2) - states(i,1) * ...
                 states(i,nsect+2) ) / objects(i,2);
      states(i,nsect+2) = states(i,nsect+2) + ...
                 dcdt(1) * delt;
%
    elseif objects(i,1) == 3,
%
% plug-flow tank
%
      m = 10 + 40 * ( 1 - objects(i,3) );
      dcdt(1) = ( states(i,2) - states(i,1) * ...
                 states(i,nsect+3-m) ) * m / objects(i,2);
      invtau = states(i,1) * m / objects(i,2);
      for j = 2:m,
        k = j - m + nsect + 2; 
        dcdt(j) = invtau * ( states(i,k-1) - states(i,k) );
      end;
      for j = 1:m,
        k = j - m + nsect + 2; 
        states(i,k) = states(i,k) + dcdt(j) * delt;
      end;
%
    elseif objects(i,1) == 4,
%
% mixing tee
%
      states(i,nsect+2) = states(i,2) / states(i,1);
%
    end;
    results(n,i+1) = states(i,nsect+2);
  end;
%
  results(n,1) = time;
  n = n + 1; 
end;
%
return;
