%  TRANSIENT RESPONSE OF NONISOTHERMAL CSTR
%  TITLE: RECDYN.M
%  DATE: November 1994
%  VERSION: 02
%  COPYRIGHT 1993, 1994 M.Gretzinger & T.Marlin
%            McMaster University   
%  This program calculates the transient response of a nonisothermal cstr
%       with a single reaction and cooling coil. It can be run without 
%       control (kc1 = kc2 = 0) or with one or two controllers. Controller 1 
%       controls the reactor temperature with the coolant flowrate. 
%       Controller 2 controls the reactor concentration with the inlet 
%       concentration. The transient response is displayed in graphical form.

%       (see Marlin's text, example 3.10 and Appendix C for derivations)

%  RECDYN.M is called by CSTR.M, the menu-driven driver program. CSTR.M 
%       has menus to load the process data (in this case contained in
%       CSTRDATA.M), modify tuning and simulation parameters and call
%       RECDYN.M.

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

clear Y1 Y2 Y3 Y4 T

%   Initialize variables
tin=tininit ;
tn=tinit ;
tcin = tcininit ;
can=cainit ;
fcbegin = fc ;    % save beginnng valve before effect of controller
cain = caininit ;

sp1 = sp1init ;
sp2 = sp2init ;
mvn=0. ;     
er1nm1 = 0 ;
er2nm1 = 0 ;

nsave = maxsave + 1 ;
time = 0. ;
cmax=round(tend/t)  ;
kont=0;

Y1=zeros(1,cmax/maxsave );
Y2=zeros(1,cmax/maxsave);
Y3=zeros(1,cmax/maxsave );
Y4=zeros(1,cmax/maxsave);
T=zeros(1,cmax/maxsave);
 
%  While loop to calculate transient response from time = 0 to time = tend.
while time < tend
    time = time + t ;
    nsave= nsave + 1;

    if time >= tdist 
         tin=tininit+tindist ;
         sp1 = sp1new ;
    end 

%  Process simulation equations **************

    can = can + t*( (f/v)*(cain-can) - ko*exp(-er/tn)*can ) ;

    tn  = tn  + t*( (f/v)*(tin-tn  ) + (a*fc^(b+1)/ ...
           (fc+ a*fc^b/(2*rhoc*cpc))) *  ...
            (1/(v*rho*cp)) * (tcin-tn) ...
           + dhrxn * ko*exp(-er/tn)*can/(rho*cp))  ;

%  Controller equations  ********************

    %  Controller 1: controls reactor temperature with coolant flowrate
    errorn1 = (sp1 - tn) ;
    if time==t, er1nm1=errorn1; end ;
    deltafc = kc1 * ( ( errorn1 - er1nm1) + ( t * errorn1)/ti1 );
    fc=fc + deltafc  ;
    if fc >= fcmax
       fc=fcmax;
    end
    if fc <= fcmin
       fc=fcmin ;   % fcmin greater than zero to prevent divide by zero
    end
    er1nm1 = errorn1 ;
    %disp(fc);

    %  Controller 2: controls reactor concentration with inlet concentration
    errorn2 = (sp2 - can);
    if time == t, er2nm1 = errorn2; end
    deltcain = kc2 * ( ( errorn2 - er2nm1) + ( t * errorn2)/ti2);
    cain = cain + deltcain;
    if cain >= 4.0
       cain=4.0;
    end
    if cain <= 0.
       cain=0. ;
    end
    er2nm1 = errorn2 ;

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

    %   Save the output variables every maxsave times to ensure student
    %   version limit on matrix size (max 1024 elements) is not exceeded.
    if nsave >= maxsave     % skip maxsave times
        kont = kont + 1;
        Y1(1,kont) = can ;  
        Y2(1,kont) = tn ;
        Y3(1,kont) = fc ;
        Y4(1,kont) = cain ;
        T(1,kont) = time ;
        nsave = 0. ;
    end

end     %while time

%  Plot results

axis('auto')
clf              % clear old graph

subplot (2,2,1), plot (T,Y1)
xlabel ('time')
ylabel ('reactor concentration')
title (['SPold = ',num2str(sp1init),', SPnew = ',num2str(sp1new),', Kc = ',num2str(kc1),', TI = ',num2str(ti1)])

subplot (2,2,2), plot (T,Y2)
xlabel('time')
ylabel('reactor temperature')

subplot (2,2,3), plot (T,Y3)
xlabel ('time')
ylabel ('coolant flowrate')
subplot (2,2,4), plot (T,Y4)
xlabel('time')
ylabel('inlet concentration')

figure (1)
pause

close all
fc=fcbegin ;   % reset fc when used in other programs