%{
Lecture 7
%}

%{
A <-> B -> C
Solve symbolically
%}

clear
home

syms t CA CA0 CB CC k1 k2 km1
[CA,CB]=dsolve('DCA=-k1*CA+km1*CB','DCB=k1*CA-km1*CB-k2*CB','CA(0)=CA0','CB(0)=0');
CC=CA0-CA-CB;

CA=simplify(CA)
CB=simplify(CB)
CC=simplify(CC)

CA0=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pseudo-steady-state case
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k1=1;
km1=1;
k2=1;

% pre-equilibrium model
CApre = @(t) CA0*exp(-(k2*k1/km1)*t);
CCpre = @(t) CA0*(1-exp(-(k2*k1/km1)*t));
% PSSA model
CApss = @(t) CA0*exp(-(k1*k2/(km1+k2))*t);
CCpss = @(t) CA0*(1-exp(-(k1*k2/(km1+k2))*t));
% full model
CAfull= @(t) subs(CA);
CCfull = @(t) subs(CC);

range=[0 20 0 1];
h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
subplot(2,1,1)
fplot(CAfull,range)
xlabel('time')
ylabel('C_A')
title('PSSA case, k_1 \approx k_{-1} \approx k_2')
hold
fplot(CApre,range,'g+')
fplot(CApss,range,'ro')
legend('Full model','Pre-equilibrium','Pseudo-steady-state')

subplot(2,1,2)
fplot(CCfull,range)
xlabel('time')
ylabel('C_C')
hold
fplot(CCpre,range,'g+')
fplot(CCpss,range,'ro')
legend('Full model','Pre-equilibrium','Pseudo-steady-state','Location','SE')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pre-equilibrium case
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k1=1;
km1=3;
k2=0.2;

% pre-equilibrium model
CApre = @(t) CA0*exp(-(k2*k1/km1)*t);
CCpre = @(t) CA0*(1-exp(-(k2*k1/km1)*t));
% PSSA model
CApss = @(t) CA0*exp(-(k1*k2/(km1+k2))*t);
CCpss = @(t) CA0*(1-exp(-(k1*k2/(km1+k2))*t));
% full model
CAfull= @(t) subs(CA);
CCfull = @(t) subs(CC);

range=[0 50 0 1];
h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
subplot(2,1,1)
fplot(CAfull,range)
xlabel('time')
ylabel('C_A')
title('Pre-equilibrium case, k_2< k_1 < k_{-1}')
hold
fplot(CApre,range,'g+')
fplot(CApss,range,'ro')
legend('Full model','Pre-equilibrium','Pseudo-steady-state')

subplot(2,1,2)
fplot(CCfull,range)
xlabel('time')
ylabel('C_C')
hold
fplot(CCpre,range,'g+')
fplot(CCpss,range,'ro')
legend('Full model','Pre-equilibrium','Pseudo-steady-state','Location','SE')
CA =
-1/2*CA0*(exp(-1/2*(k1+km1+k2-(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*k1-exp(-1/2*(k1+km1+k2-(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*km1-exp(-1/2*(k1+km1+k2-(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*k2-exp(-1/2*(k1+km1+k2-(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2)-exp(-1/2*(k1+km1+k2+(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*k1+exp(-1/2*(k1+km1+k2+(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*km1+exp(-1/2*(k1+km1+k2+(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*k2-exp(-1/2*(k1+km1+k2+(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))/(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2)
CB =
CA0*k1*(exp(-1/2*(k1+km1+k2-(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)-exp(-1/2*(k1+km1+k2+(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t))/(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2)
CC =
-1/2*CA0*(-2*(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2)+exp(-1/2*(k1+km1+k2-(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*km1+exp(-1/2*(k1+km1+k2-(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2)-exp(-1/2*(k1+km1+k2+(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*km1+exp(-1/2*(k1+km1+k2+(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2)+exp(-1/2*(k1+km1+k2-(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*k1+exp(-1/2*(k1+km1+k2-(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*k2-exp(-1/2*(k1+km1+k2+(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*k1-exp(-1/2*(k1+km1+k2+(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2))*t)*k2)/(k1^2+2*k1*km1-2*k1*k2+km1^2+2*k2*km1+k2^2)^(1/2)
Current plot held
Current plot held
Current plot held
Current plot held