%============================================================
%             |  Author: James C. Sutherland  |
%             |            CHEN 1703          |
%             |-------------------------------|
%
% This script performs regression on data to determine the
% rate constant for the reaction discussed in class.
% See the lecture notes for more information on the
% formulation of this problem.
%
%============================================================
 
clear; clc; close all;

% define the temperature and rate constant data.
T = [ 313 319 323 328 333 ];
k = [ 0.00043 0.00103 0.00180 0.00355 0.00717 ];

m = length(T);
R = 8.314;  % gas constant, J/mol-K

%--------------------------------- 
% METHOD 1: Use MATLAB's polyfit to do the regression.
% Here we are fitting to a first-order (linear) polynomial.
aa = polyfit( -1./(R*T), log(k), 1 );
fprintf('From polyfit: A=%1.2e, Ea=%1.2e\n',exp(aa(2)),aa(1));
%--------------------------------- 

%--------------------------------- 
% METHOD 2:
A = [ ones(m,1) -1./(R*T') ];
b = log(k');

% form and solve the normal equations
AA = A'*A;
bb = A'*b;
xx = AA\bb;

fprintf('From normal equations: A=%1.2e, Ea=%1.2e\n',exp(xx(1)),xx(2));
%--------------------------------- 


% plot some stuff
tmp = -1./(R*linspace(312,334));
ytmp = xx(2)*tmp + xx(1);
plot(-1./(R*T),log(k),'ko',tmp,ytmp,'r-');
xlabel('-1/RT (mol/J)');
ylabel('log(k) (1/s)');
legend('data','best fit');
text(-3.85e-4,-5,strcat('ln(A)=',num2str(xx(1))));
text(-3.85e-4,-5.25,strcat('Ea=',num2str(xx(2))));
