%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Batch Reactor Data Analysis (Zero Order)
% Warren Stiver, January 2016
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This script fits Concentration vs time data to a Zero Order reaction
% kinetic model.
% The resulting graphical representation includes the best fit to the data, 
% and the uppper / lower limits of the fit based on a specified confidence limit
%
% The best estimate of the reaction rate coefficient & the upper and lower
% estimate are determined. 
%
% Statistics reference:  
% Kirkpatrick E.G. (1974) Introductory Statistics and Probability for 
%       Engineering, Science, and Technology, Prentice-Hall, 446p.
%       Chapter 10.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  DATA

%C = [100 95 90 85 80 75 70 65 60 55 50]   % straight line for checking!
C = [100 92 91 85 80 78 72 63 60 54 50];    % g/mL
t = [0 1 2 3 4 5 6 7 8 9 10];               % s


%%%%%%%%%%%%%%%%%%%%%
% desired confidence limits

CLdesired = 0.95; % 95% confidence limits for uncertainty of fitted line


%%%% regression using Matlab's PolyFit function
zeroOrder = polyfit(t, C, 1);   % C as a function of t with 1 corresponding to the degree of the polynomial (highest power on t)
C_fit = zeroOrder(1)*t+zeroOrder(2); % (1) is slope with (2) as the intercept

%%%% regression from fundamental relationships
n = length(t);
t_avg = mean(t);
C_avg = mean(C);
numerator = 0;
denominator = 0;

for i= 1:n
    numerator = numerator + ((t(i)-t_avg) * (C(i)-C_avg));
    denominator = denominator + (t(i)-t_avg)*(t(i)-t_avg);
end
slope = numerator/denominator;
intercept = C_avg - slope * t_avg;

C_man = slope * t + intercept;  %man ... manually fit

%%%% determining the regression coefficient
devRnum = 0;
devRden = 0;
for i= 1:n
    devRnum = devRnum + power(C(i)-C_man(i),2);
    devRden = devRden + power(C(i)-C_avg, 2);
end

Rsquared = 1-devRnum/devRden   % the degree to which the fit explains the difference between C measured and the avg C

%%%% determining the uncertainty in the regression
devY = 0;
devX2 = 0;
sumX2 = 0;
for i = 1:n
    devY = devY + power(C(i)-C_man(i),2);
    devX2 = devX2 + (t(i)-t_avg)*(t(i)-t_avg);
    sumX2 = sumX2 + t(i)*t(i); 
end

Syx = power( devY /(n-2),0.5);              % regression uncertainty
Sb = Syx/(power(devX2, 0.5));               % slope uncertainty
Sa = Syx * (power(sumX2/n/devX2, 0.5));     % intercept uncertainty

tstat = tinv((1+CLdesired)/2,n-2);          
        % tinv is an embedded MatLab function to determine the t
        % distribution value for a given statistical probability and
        % degrees of freedom (number of data points minus 2)

slope_upper = slope + tstat*Sb;
slope_lower = slope - tstat*Sb;

k = slope                          % zero order units of g mL^-1 s^-1 (or g/mL/s) based on C and t units
kupper = slope_upper
klower = slope_lower


for i = 1:n
    C_upper(i) = intercept + slope*t(i) + tstat * (Syx * power((1/n + power((t(i)-t_avg),2)/devX2),0.5));
    C_lower(i) = intercept + slope*t(i) - tstat * (Syx * power((1/n + power((t(i)-t_avg),2)/devX2),0.5));
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Plotting data in various forms
%
% figure with Measured data (symbol only) and Fitted line (line only)
figure(1)
plot (t,C, 'b*')
hold on
grid on
plot(t, C_fit, 'b')

xlabel('Time (s)', 'FontSize', 16);
ylabel('Concentration (g/mL)', 'FontSize', 16);
hleg = legend('Measured', 'Fit MatLAB');
fh = figure(1);
set (fh, 'Color', 'w');
%axis([0 2500 0 0.002])


%%%%
% figure with Measured data (symbol only) and Fitted line using base equations (line only)
figure(2)
plot (t,C, 'k*')
hold on
plot(t, C_man, 'k')

xlabel('Time (s)', 'FontSize', 16)
ylabel('Concentration (g/mL)', 'FontSize', 16)
hleg = legend('Measured',  'Manual Fit');
str = ['C = ', num2str(slope), ' t + ', num2str(intercept)];
annotation('textbox', [0.55 0.7 0.3 0.05],'String', str, 'Fontsize', 16, 'FitBoxtoText', 'on');
fh = figure(2);
set (fh, 'Color', 'w');
%axis([0 2500 0 0.002])

%%%%
% figure with Measured data (symbol only), Best, Upper, Lower fit lines 
figure(3)
plot (t,C, 'rd', 'MarkerFaceColor', 'r', 'markers', 12 )
hold on

plot(t, C_man, 'r')
plot(t, C_upper, 'r--');
plot(t, C_lower, 'r--');

xlabel('Time (s)', 'FontSize', 16)
ylabel('Concentration (g/mL)', 'FontSize', 16)
hleg = legend('Measured',  'Best Fit', 'Upper/Lower Limits');
fh = figure(3);
set (fh, 'Color', 'w');
%axis([0 2500 0 0.002])
