%  *****  EXAMPLE 3.3  SERIES CSTR's   *******
%  *****  from MARLIN
%  COPYRIGHT 1994 by McMASTER UNIVERSITY
%  WRITTEN BY T. MARLIN
%  FILE NAME : EXSINE.M (EXAM3_3.M)
%  LAST UPDATE: NOVEMBER 1994
%  VERSION 1.0
%  DESCRIPTION OF FILE:
%  this m-file calculates the transient response 
%  for sine input forcing to CSTR's from example 3_3
%  *********************************************

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

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

%  *********************************************
%   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   ;
omega1 = .1   ;  % initial frequency of sine forcing
omega2 = .5   ;  % modified frequency of sine forcing
amp  = 0.5    ;  % amplitude of sine forcing
%  ************************************
%      SIMULATION DATA
%  ************************************
delt = 0.5   ;  % step size
tstart=0.    ;
tend=900.    ;
t1 = 300     ;  % time of first change in frequency
t2 = 600     ;  % time of return to original frequency

cmax=round(tend/delt)  ;
nskip = 2   ;  % 1/nskip is fraction of results saved for plot
cntskip = 0 ;  % initialize
nstore = 1  ;  % initialize
%  ************************************
%    INITIAL CONDITIONS AND STORAGE VECTORS
%  *************************************
cain = .925        ;
ca1init=0.616     ; 
ca2init = 0.41   ;
ca1=ca1init     ;
ca2=ca2init   ;
%
Y=zeros(1,cmax/nskip);
MV=zeros(1,cmax/nskip);
T=zeros(1,cmax/nskip);
Y2=zeros(1,cmax/nskip) ;
% **************************************
%  SIMULATION
% **************************************

for cnt=1:1:cmax
%
%  introduce the sine forcing in concentration

   ca0 = cain + amp* sin (omega1*cnt*delt);
   if cnt*delt > t1
        ca0 = cain +  amp * sin (omega2*cnt*delt);
   end
   if cnt*delt > t2
        ca0 = cain + amp * sin (omega1*cnt*delt) ;
   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;
%
    cntskip = cntskip + 1 ;
    if cntskip >= nskip
       Y(1,nstore) = ca1      ;
       Y2(1,nstore) = ca2     ;
       MV(1,nstore)=ca0       ;
       T(1,nstore) = cnt*delt ;
       YLIMIT(1,nstore) = .85 ;
       cntskip = 0.;
       nstore = nstore + 1 ;
    end
%
end
%
%  ******************************************
%     PLOT DATA
%  ******************************************
subplot (2,2,1), plot (T(1,:),Y(1,:))
xlabel ('time')
ylabel ('reactor 1 concentration')
title ('time-domain plot of frequency response')
%
subplot (2,2,3), stairs (T(1,:),MV(1,:))
xlabel('time')
ylabel('feed concentration')
%
subplot (2,2,2), plot (T(1,:),Y2(1,:))
ylabel ('reactor 2 concentration')

figure (1) ; pause ; close all