%% ... Rate constant estimation
%% ... Reaction:
%% ... CO + Cl2 --> COCl2
clc;
%% ... Experimental data obtained in constant volume BATCH
% ... Time [s]
time = [0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38]';
% ... Pressure p(t) [Pa]
p = [599554	434817 393038 373953 378726 357128 351020 336651 ...
    343442	337214	325668	333051	342957	327858	321577	325457 ...
    334276	319888	329878	317612]';
[ndata,null] = size(time);
XC = zeros(ndata);
% ... Initial CO molar fraction
yA0 = 0.5;
p0  = 599554; % [Pa]
% ... Other values
R   = 8.314;         % [J/mol/K]
T   = 423;           % [K]
cA0 = p0*yA0./(R*T); % [mol/m3]
%% ... Conversion and integral
XA = (p0 - p)./(yA0.*p0);
IX = @(X)(2.*sqrt(1-X))./(3.*X.^2-6.*X+3)-2./3;


%% ... Solving the system t^T*IX = (t*t^T)*K
TT = (time'*time); 
%
K = TT\(time'*IX(XA));
% ... Estimate error by computing the covariance "matrix"
err = sqrt(TT^(-1));
%[K,stdx,mse] = lscov(TT,time'*IX);
k_alpha = K ./( cA0.^(3/2));
fprintf('Rate constant for catalyst alpha\n')
fprintf('k_alpha = %10.4e +- %10.4e (m3/mol)^3/2 /s\n',k_alpha,err);
t = 0;
FCN = @(X,t)((2.*sqrt(1-X))./(3.*X.^2-6.*X+3)-2/3 - K .*t);
fprintf('\n%12s %12s %12s\n','time','Xcalc','Xexpt');
for i = 1 : ndata
    t     = time(i);
    XC(i) = fzero(@(x)FCN(x,t),XA(i));
    fprintf('%12.2f %12.4f %12.4f\n',t, XC(i), XA(i));
end
subplot(2,1,1)
plot(time,IX(XA),'*',time,K.*time);
xlabel('Time / s')
ylabel('Integral () ')
legend('Data','Fit')

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