%  *****  EXAMPLE 3.3  SERIES CSTR's   *******
%  *****  from MARLIN
%  COPYRIGHT 1994 by McMASTER UNIVERSITY
%  WRITTEN BY T. MARLIN
%  FILE NAME : EXAM3_3.M
%  LAST UPDATE: JULY 1994
%  VERSION 1.0
%  DESCRIPTION OF FILE:
%  this m-file calculates the transient response 
%  for step to CSTR, examle 3_3
%  *********************************************

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

%  *********************************************

%  data 
%
%  *********************************************
%   PROCESS DATA
%  *********************************************
v1 = 1.05     ; % volume of tank 1
v2 = 1.05     ; % volume of tank 2
k = 0.040     ; % reaction rate constant
f = 0.085     ; % flow rate
ca0 = 0.925   ;

tau=8.23        ; % process time constant (both)
kp=0.669       ; % process gain (both)
%beta = kp*(1-exp(-t/tau)) ;
%alpha = exp(-t/tau)       ;
%  ************************************
%      SIMULATION DATA
%  ************************************
delt = 0.1     ;  % step size
tstart=0.    ;
tend=60.    ;
%  ************************************
%    INITIAL CONDITIONS AND STORAGE VECTORS
%  *************************************
ca1init=0.616     ; 
ca2init = 0.41   ;
ca1=ca1init     ;
ca2=ca2init   ;

cmax=round(tend/delt)  ;
%
Y=zeros(1,cmax);
MV=zeros(1,cmax);
T=zeros(1,cmax);
Y2=zeros(1,cmax) ;
YLIMIT=zeros(1,cmax) ;

% **************************************
%  SIMULATION
% **************************************

for cnt=1:1:cmax
%
% ******  LINEARIZED EXPRESSION  ******    

    YANAL(cnt) = ca2init;
    ca0 = 0.9250 ;
    if cnt*delt > 10
          ca0 = 2*0.925 ;
          YANAL(cnt) = 0.828 - 0.050*(cnt*delt-10)...
                      *exp (-(delt*cnt-10)/8.25)...
                  +(0.414-0.828)*exp(-(delt*cnt-10)/8.25) ;
    end

% ******  NON-LINEAR NUMERICAL SIMULATION  *********

   ca1dot =  (f*(ca0 - ca1) - v1*k*ca1)/v1 ;
   ca2dot =  (f*(ca1 - ca2) - v2*k*ca2)/v2 ;
   ca1 = ca1 + delt * ca1dot;
   ca2 = ca2 + delt * ca2dot;
%
    Y(1,cnt) = ca1      ;
    Y2(1,cnt) = ca2     ;
    MV(1,cnt)=ca0      ;
    T(1,cnt) = cnt*delt   ;
    YLIMIT(1,cnt) = .85;
%
end
%
%  ******************************************
%     PLOT DATA
%  ******************************************

subplot (2,2,1), plot (T(1,:),Y(1,:))
xlabel ('time')
ylabel ('tank 1 concentration')
axis([0 tend .30 1.3 ])
title ('Example 3.3')
%
%
subplot (2,2,3), stairs (T(1,:),MV(1,:))
axis ([0 tend .5 2 ])
xlabel('time')
ylabel('inlet concentration')
%
subplot (2,2,2), plot (T(1,:),Y2(1,:),T(1,:),YLIMIT(1,:),'--',...
                      T(1,:), YANAL(:),':')
axis ([0 tend .3 1.3])
ylabel ('tank 2 concentration')

figure (1)
pause
close all