%% ... Rate constant estimation
%% ... Reaction:
%% ... CO + Cl2 --> COCl2
clc;
%% ... Experimental data obtained in isothermal, isobaric PFR
% ... Residence time [s]
tr = [45.9 22.9 9.2 4.6 2.3 1.5 1.1 0.9 0.8 0.7 0.5 0.2]';
% ... Conversion
XA = [0.997 0.995 0.992 0.987 0.979 0.971 0.965 0.958 0.952 0.946 0.928 0.871]';
[ndata,null] = size(tr);
IXA = zeros(ndata,1);
XC  = zeros(ndata,1);
eps = zeros(ndata,1);
info = zeros(ndata,1);
% ... Initial CO molar fraction
yA0 = 0.5;
p   = 300000; % [Pa]
% ... Other values
R   = 8.314;         % [J/mol/K]
T   = 423;           % [K]
cA0 = p.*yA0./(R*T); % [mol/m3]
VR  = 2e-6;          % [m3]

%% ... Conversion and integral
IX = @(X)(((1-yA0.*X)./(1-X)).^(5/2));

for i = 1 : ndata
    IXA(i) = quadl(IX,0,XA(i));
end

%% ... Solving the system t^T*IX = (t*t^T)*K
TT = (tr'*tr); 
% ... Matrix approach
K = TT\(tr'*IXA);
% ... or using std matlab fcn
%[K,stdx,mse] = lscov(TT,tr'*IXA);
% ... Estimate error by computing the covariance "matrix"
err = sqrt(TT^(-1));
k_beta = K ./( cA0.^(3/2));
fprintf('Rate constant for catalyst beta\n')
fprintf('k_beta = %10.4e +- %10.4e (m3/mol)^3/2 /s\n',k_beta,err);
t = 0;
FCN = @(X,t)(( IX(X) - K .*t) );
tol = 1e-9;
fprintf('\n%12s %12s %12s %12s %4s\n','tr','Xcalc','Xexpt','eps','Iter');
for i = 1 : ndata
    t     = tr(i);
    %[XC(i),eps(i),info(i)] = newrap(@(x)FCN(x,t),0.99*XA(i),tol);
    
    [XC(i),eps(i),info(i)] = fsolve(@(x)FCN(x,t),0.99*XA(i));
    
    fprintf('%12.2f %12.4f %12.4f %12.4e %4i\n',t, XC(i), XA(i),eps(i), info(i));
end

subplot(2,1,1)
plot(tr,IXA,'*',tr,K.*tr);
xlabel('Residence time / s')
ylabel('Integral () ')
legend('Data','Fit')

subplot(2,1,2)
plot(tr,XA,'*',tr,XC);
xlabel('Time / s')
ylabel('Conversion')
legend('Data','Fit')

