% *****  EXAMPLE 5.4  SERIES OF TWO CSTR's, INVERSE RESPONSE
% *****             from  MARLIN               ************
% COPYRIGHT: 1994 by McMASTER UNIVERSITY
% WRITTEN BY : T. MARLIN
% FILENAME: EXAM5_4.M
% LAST UPDATED: November 1994
% VERSION 2.0
% DESCRIPTION OF FILE:
%  this m-file calculates the transient response 
%  for example 5_4 with numerical solution to ode
%  and analytical solution
%  **********************************************************

clear ; clf ; axis ('auto') ; clc

% **************************************************************
%   PROCESS DATA
% *************************************************************

nstage = 2 ;         % number of reactors
v = 1.05        ;
f = 0.085        ;
k = 0.040      ;
ca(1) = 0.925    ;   % feed concentration
 
for nn=1:nstage      %  initialize the tank concentrations
  ca(nn+1) = ca(1)*(f/(f+v*k))^nn ;
end
%  ***************************************************
%  SIMULATION DATA & SIZING ARRAYS
%  ***************************************************
delt=.2         ;
tstart=0.    ;
tend=60.    ;
%
mvn=0       ;
cmax=round(tend/delt)  ;
%
Y=zeros(1,cmax);
MV=zeros(1,cmax);
T=zeros(1,cmax);
Y1=zeros(1,cmax);
Y2=zeros(1,cmax) ;
% *****************************************************
%  SIMULATION CALCULATIONS
% *****************************************************

for cnt=1:1:cmax
%                          enter the step input and calculate
    canal(cnt) = 0. ;      % the linearized analytical solution
    if cnt*delt > 10
          f   = .085 + .0085 ;     % step change
          ca(1) = 0.925 - .085 ;
                                   % analytical step response
          canal(cnt) = -0.0141 +( .0141 +.00337*(cnt*delt-10)) ...
                         *exp (-(cnt*delt-10)/8.25) ;
    end
    
    for nn=1:nstage       % simulate the non-linear model
                          % by integrating numerically
     ca(nn+1) = ca(nn+1) + delt *( f*(ca(nn)-ca(nn+1))/v - k*ca(nn+1))  ;
    end
% ******************************************************************
%  STORE RESULTS IN ARRAYS
% *******************************************************************
    Y(1,cnt) =  ca(2)     ;
    Y2(1,cnt) = ca(3)    ;
    MV(1,cnt)=  f      ;
    T(1,cnt) = cnt*delt   ;
%
end
% ******************************************************************
% PLOT RESULTS
% ******************************************************************

subplot (2,1,1), plot (T(1,:),Y2(1,:), T(1,:),canal(:)+0.4145,'--')
ylabel ('tank 2 concentration')
axis ([0 tend .3900 0.43])
title('Example 5.4, 2 Series reactors (dashed is linearized approx.)')

subplot (2,1,2), stairs (T(1,:),MV(1,:))
axis ([0 tend .08 .1 ])
xlabel('time')
ylabel('solvent flow')

%

figure (1)
pause

close all