% 1.731 Chemical reactor example
function cstr
% Program to simulate response of cstr to variable inflow solute concentration
% Conc units are mgl^(-1), flow units are m^3/hr,
% volume units are m^3, decay coef. units are hr^(-1)
% Load inflow solute concentration data
load inconc.dat
n=length(inconc)
%define inputs
% specify time step (hrs)
h=100
%specify flow rate (m^3/hr), volume (m^3), initial conc (g/m^3)
%decay rate(in hr^(-1))
Q=1.
V=1000.
conc(1)=0.
amax=0.001
concmax=4
% time loop (increment t from 1 through n)
for t=1:n
   % compute nominal storage (in volume units) at time t+1,
   % for given storage at time t
   % derive nonlinear decay rate from most recent concentration
   a=-amax*conc(t)*(conc(t)-concmax)/concmax;
   a=a*((conc(t)>0)&(conc(t)<concmax));
   conc(t+1)=conc(t)+(Q*h/V)*inconc(t)-(Q*h/V)*conc(t)-a*h*conc(t);
end
% print results
inconc
conc
% close all old plot windows and set up new plot window
close all
figure
% construct conc plot
subplot(2,1,1)
bar(h*[0.5:n-0.5],inconc)
subplot(2,1,2)
plot(h*[0:n],conc)
return