%{
 Lecture 8, Analyzing Reactor Kinetics Data
%}

%
%  m-xylene bromination
%
clear;
home;

t=[0 2.25 4.50 6.33 8 10.25 12 13.5 15.6 17.85 19.6 27 30 38 41 45 47 57 63]'; % min
CBr2=[.3335 .2965 .2660 .2450 .2255 .2050 .1910 .1794 .1632 .1500 ...
    .1429 .116 .1053 .0830 .0767 .0705 .0678 .0553 .0482]'; % M

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
subplot(2,1,1);
plot(t,CBr2,'ro'); hold;
p=polyfit(t,CBr2,4);  % fit 4th order polynomial to data
tp=[0:1:65];
f=polyval(p,tp);       % evaluate polynomial
plot(tp,f);
legend('C_{Br_2}','Polynomial fit')
title('Xylene catalytic bromination experimental data')
xlabel('Time (min)')
ylabel('Concentration (M)')

%
%  rate = k C^alpha
%  log(rate)=log(k)+alpha*log(C)
%
k=polyder(p);          % differentiate polynomial C(t)
fp=polyval(k,t);     % evaluate derivative C'(t)
lfp=log(-fp);          % take log(-C'(t))
lCBr2=log(CBr2);
subplot(2,1,2);
plot(lCBr2,lfp,'r^')
xlabel('ln(C_{Br_2})')
ylabel('ln(-dC_{Br_2}/dt)')
%
%  Note incorrect curvature at one endpoint
%

% need to remove last data point
lCBr2=lCBr2(1:18); lfp=lfp(1:18);

X=[ones(size(lCBr2)) lCBr2];
A=X\lfp;  % this does linear regression of data against log(CBr2)
R=corrcoef(lCBr2,lfp);

k=exp(A(1));
alpha=A(2);

sprintf('Reaction order in Br2=%f, rate constant=%f, correlation coef=%f',alpha,k,R(1,2))
%fit=polyval(A,lCBr2);
fit=A(1)+A(2)*lCBr2;
hold;
plot(lCBr2,fit)
legend('Rates from polynomial differentiation','Linear fit','Location','SE')

%
% fit data to integrated rate law
%
%  (1/2)kt=(CBr2^-0.5-CBr20^-0.5)
%
sqrtCBr2=CBr2.^(-0.5);
h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
plot(t,sqrtCBr2,'ro')
xlabel('Time (min)')
ylabel('C_{Br_2}^{-1/2}')
X=[ones(size(t)) t];
A=X\sqrtCBr2;
R=corrcoef(t,sqrtCBr2);
fit=A(1)+A(2)*t;
hold;
plot(t,fit);
legend('Concentration data','Linear fit','Location','SE')

k=2*A(2);
b=1/A(1)^2;
sprintf('Rate constant=%f/s M^(1/2), intercept=%f, correlation coefficient=%f',k,b,R(1,2))
title('Fit of xylene bromination data to integrated rate law')
Current plot held
Warning: Imaginary parts of complex X and/or Y arguments ignored.
ans =
Reaction order in Br2=1.521858, rate constant=0.095362, correlation coef=0.996237
Current plot held
Current plot held
ans =
Rate constant=0.089160/s M^(1/2), intercept=0.325648, correlation coefficient=0.999726