%{
  Lecture 14: Non-isothermal CSTR steady states
%}
clear;
home;

% adiabatic CSTR reaction, A -> B
R0 = 8.31441;

Ea = 125000;  % J/mol
k0 = 2.9e20;  % min-1
dH = -90000;  % J/mol
rhoCp = 4000; % J/l K
CA0 = 2;      % M
tau = 1;      % min
kappa = 0;
J = -(CA0/rhoCp)*dH;
Ta = 273;

k= @(T) k0*exp(-Ea./(R0*T));
X= @(T) tau*k(T)./(1+tau*k(T));

%T=[280:1:350];
%G=X(T);
limits=[280 350];

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
hold all % prevents color from being reset to blue each time
fplot(X,limits);
xlabel('T (K)'); ylabel('R(T) and G(T)')
axis([280 350 0 1])
title({'A \rightarrow B CSTR Multiple Steady States';'T_0 = 290 - 306 K'})

i=figure('color',[1 1 1]);
set(i,'WindowStyle','Docked')
hold

m=(1+kappa)/J;
for T0=290:1:306   % create the R lines
    b=-(kappa*Ta+T0)/J;
    R=@(T) m*T+b;
    figure(h)
    if mod(T0,2)==0; fplot(R,limits); end

    f=@(T) X(T)-R(T);
    for T1=280:10:350  % search for intersections with G
        TS=fzero(f,T1);
        if mod(T0,2)==0; figure(h); plot(TS,R(TS),'ro'); end
        figure(i); plot(T0,TS,'ro')
    end
end
figure(i);
xlabel('T_0 (K)'); ylabel('T_{ss} (K)'); axis([290 306 290 350])
title({'A \rightarrow B CSTR Multiple Steady States';'Steady state Temperatures vs. T_0'})

% Set T0=300, vary CA0, G(T)=X(T) is constant, R(T) varies
T0=300;
j=figure('color',[1 1 1]);
set(j,'WindowStyle','Docked')
hold all % prevents color from being reset to blue each time
limits=[295 370];
fplot(X,limits);
xlabel('T (K)'); ylabel('R(T) and G(T)')

for CA0=1:.25:3
    J = -(CA0/rhoCp)*dH;
    m=(1+kappa)/J;
    b=-(kappa*Ta+T0)/J;
    R=@(T) m*T+b;

    fplot(R,limits)

    f=@(T) X(T)-R(T);
    for T1=300:10:380  % search for intersections with G
        TS=fzero(f,T1);
        plot(TS,R(TS),'ro')
    end
end
axis([295 375 0 1])
title({'A \rightarrow B CSTR Multiple Steady States';'C_{A0} = 1 - 3 M'})
Current plot held