% LEP 12-6 : Series Reactions in a CSTR
%

tspan = [0 300]; % Range for the independent variable time 
y0 = [273.]; % Initial values for the dependent variable T
Cp=200;
Cao=0.3;
To=283;
DH1=-55000;
DH2=-71500;
vo=1000;
E2=27000;
E1=9900;
UA=40000;
Ta=330;
V=10;
[t, y]=ode45(@(t,y)ODEfun(t,y,Cp,Cao,To,DH1,DH2,vo,E2,E1,UA,Ta,V),tspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,600,400])
kappa = UA / (vo * Cao) / Cp;
Tc = (To + kappa * Ta) / (1 + kappa); 
tau = V/vo; 
z=size(t);
for i=1:z(1,1)
    k2(i) = 4.58 * exp(E2 / 1.987 * (1 / 500 - (1 / y(i,1)))); 
    k1(i) = 3.3 * exp(E1 / 1.987 * (1 / 300 - (1 /y(i,1)))); 
    R(i)=Cp*(1+kappa)*(y(i,1)-Tc);
    G(i)=(-DH1*tau*k1(i)/(1+tau*k1(i)))-(tau*k1(i)*tau*k2(i)*DH2)/(1+tau*k1(i))/(1+tau*k2(i));
    Ca(i) = Cao / (1 + tau * k1(i)); 
    Cb(i) = tau * k1(i) * Ca(i) / (1 + k2(i) * tau); 
    rB(i)=k1(i)*Ca(i)-k2(i)*Cb(i);
    Cc(i)=Cao-Ca(i)-Cb(i);
    Sinsta(i)=rB(i)/(k2(i)*Cb(i));
    Soverall(i)=Cb(i)/(Cc(i)+0.0000001);
end
nexttile
plot(y(:,1),G(1:z(1,1))/100,y(:,1),R(1:z(1,1))/100,'LineWidth',1.5);
leg=legend('G(T)','R(T)');
ylabel('G(T),R(T) *10^2 (J/mol)');
xlabel('T(K)');
axis([250 800 0 1500])
title('Heat curves')
leg.ItemTokenSize=[10,10];
grid on
nexttile
plot(y(:,1),Cb(1:z(1,1))*1000,'LineWidth',1.5);
leg=legend('C_B');
ylabel('C_B*10^-3(mol/dm^3)');
xlabel('T(K)');
axis([250 800 0 600]);
leg.ItemTokenSize=[10,10];
title('Plot of C_B')
grid on
nexttile
plot(y(:,1),Soverall(1:z(1,1)),'LineWidth',1.5)
leg=legend('S_{B/C}');
leg.ItemTokenSize=[10,10];
title('Selectivity')
ylabel('Overall Selectivity');
xlabel('T(K)');
axis([250 800 0 1000]);
grid on
nexttile
plot(y(:,1),Sinsta(1:z(1,1)),'LineWidth',1.5)
leg=legend('S_{B/C}');
leg.ItemTokenSize=[10,10];
title(' Selectivity')
ylabel('Instantaneous Selectivity');
xlabel('T(K)');
axis([250 800 0 100000]);
grid on

%

function dYfuncvecdt = ODEfun(t,Yfuncvec,Cp,Cao,To,DH1,DH2,vo,E2,E1,UA,Ta,V); 
T = Yfuncvec(1); 
% Differential equation
dTdt = 2; 
dYfuncvecdt = [dTdt]; 
end