%% **** Program eaestim.m
%% ****
%% **** Linear (or linearized) regression of k(T) data
%% ... Read input data file ( T | k(T)) 
Rg = 8.314; % J/mol/K
data = dlmread('data.txt');
%% ... Transform input data (linearize it)
x    = 1./(data(:,1));
y    = log(data(:,2));
[s,t]= size(data);
stdv = ones(s,1);
if (t==3)
    stdv=data(:,3);
end
w  = 1./stdv.^2;
S  = sum(w);
Sx = sum(w.*x);
Sy = sum(w.*y);
Sxx= sum(w.*x.^2);
Sxy= sum(w.*x.*y);

Delta = S*Sxx - (Sx)^2;
a1    = (Sxx*Sy - Sx*Sxy)./Delta;
a2    = (S*Sxy - Sx*Sy)./Delta;
y_fit = a1 + a2*x;
EA    = -a2*Rg;
fprintf('\n Activation energy = %6.2f kJ/mol \n',EA/1e3)
fprintf('\n ln(kinf)             = %6.2f \n',a1)

if (q==t)
    h=plot(x,y,'rs','MarkerFaceColor','r');
    title('Unweighted fit')
else
    clf;set(gcf,'color','w');
    h=errorbar(x,y,stdv,'rs','MarkerFaceColor','r');
    title('Weighted fit')
end
hold on;
q=plot(x,y_fit,'b.--','linewidth',2);
xlabel('1/T [1/K]')
ylabel('ln(k)')
legend([h(1),q(1)],'Data',sprintf('\nEa = %6.2f kJ/mol \nln(kinf) = %6.2f \n',EA/1e3,a1),'location','Northeast')
