% LEP 10-6 : Decay in a Straight-Through Transport Reactor
%

zspan = [0 10]; % Range for the independent variable z
y0 = [0]; % Initial values for the dependent variables i.e X
Ka=0.05;
Kb=0.15;
Pao=12;
eps=1;
A=7.6;
T=673;
rho=80;
kprime=0.0014;
D=1.5;
Uo=2.5;
Kc=0.1;
[z,y]=ode45(@(z,y)ODEfun(z,y,Ka,Kb,Pao,eps,A,T,rho,kprime,D,Uo,Kc),zspan,y0);
s=size(z);
tiledlayout(1,2)
set (gcf,'Position',[0,0.1,800,300])
for i=1:s(1,1)
    U(i)=Uo*(1+eps*y(i));
    a(i)=1/(1+A*(z(i)/U(i))^0.5);
    Pa(i)=Pao*(1-y(i))/(1+eps*y(i));
    Pb(i)=Pao*y(i)/(1+eps*y(i));
    Pc(i)=Pb(i);
    raprime(i)=a(i)*(-kprime*Pa(i)./(1+Ka*Pa(i)+Kb*Pb(i)+Kc*Pc(i)));
    ra(i)=rho*raprime(i);
end
nexttile
plot(z,y,z,a(1:s(1,1)),'LineWidth',2)
ylabel('X, a');
xlabel('z');
leg=legend('X','a');
leg.ItemTokenSize=[10,10];
axis([0 10 0 1])
title('Activity and conversion profiles')
grid on
nexttile
plot(z,-ra(1:s(1,1)),'LineWidth',2)
ylabel('Rate');
xlabel('z');
leg=legend('-r_A');
leg.ItemTokenSize=[10,10];
axis([0 10 0 0.5])
title('Rate profile')
grid on

%

function dYfuncvecdZ = ODEfun(z,Y,Ka,Kb,Pao,eps,A,T,rho,kprime,D,Uo,Kc); 
X=Y(1);
% Explicit equations
R=0.082;
U=Uo*(1+eps*X);
Pa=Pao*(1-X)/(1+eps*X);
Pb=Pao*X/(1+eps*X);
Cao=Pao/R/T;
Pc=Pb;
a=1/(1+A*(z/U)^0.5);
raprime=a*(-kprime*Pa/(1+Ka*Pa+Kb*Pb+Kc*Pc));
ra=rho*raprime;
% Differential equations
dXdz = -ra/Uo/Cao;
dYfuncvecdZ = [dXdz]; 
end
