% Example of differentiating experimental data
% by fitting polynomial to groups of data and differentiating
% the polynomial
%
% HIT SPACE BAR TO STEP THROUGH PROCESS ON PLOT
 
clear
fprintf('--------------------------\n') % run separator in command window
 
% Differential method of analysis of batch reactor data in ReactorLab
% copy data to Matlab file, comment the header lines, put all data into 
% array with arbitrary name d
 
% D1L1 Batch, Data Set 2b
% 
% k*_2b remained constant at    1.00E-4  (value at 300 K)
% Ea_2b remained constant at       60.0  kJ/mol
% n_2b remained constant at       2.00  (dimensionless)  
% T_2b remained constant at      325.0  K
% Cin_2b remained constant at      100.0  mol/m3
% vol_2b remained constant at       1.00  m3
% time_2b has units of s
% tau_2b has units of s
% Cout_2b has units of mol/m3
% conv_2b has units of (dimensionless)
 
% time_2b   tau_2b  Cout_2b conv_2b
d = [0.00   0.00    100.0   0.00
1.00    1.00    95.8    4.21E-2
2.00    2.00    87.1    0.129
3.00    3.00    83.1    0.169
4.00    4.00    80.4    0.196
5.00    5.00    74.9    0.251
7.50    7.50    68.2    0.318
10.0    10.0    61.0    0.391
20.0    20.0    43.8    0.562
30.0    30.0    34.8    0.652
40.0    40.0    27.9    0.721
50.0    50.0    23.8    0.762
60.0    60.0    20.3    0.797
70.0    70.0    18.4    0.816
80.0    80.0    16.4    0.836
90.0    90.0    14.7    0.853];
 
% extract time and concentration vectors
t = d(:,1);
c = d(:,3);
 
% plot
plot(t,c,'.')
title('HIT SPACE KEY TO STEP - differentiating data')
ylabel('conc (mol/m3)')
xlabel('time (s)')
 
[rows cols] = size(c); % size is std matlab function
npts = rows;
 
%----- DIFFERENTIATE WITH GENERAL POLYNOMIAL FIT ----------
 
% We want to look at the rate of change of the concentration vs. t
% by differentiating these data.
% The simplest way is just to get the slope between each
% pair of points using Matlab's diff function.
% But this will show a very "noisy" derivative with these data.
% We can also fit a polynomial to groups of points
% to get some "smoothing" of the data.
 
% The following will also work with pts that are
% NOT equally spaced in time
 
% Vary n and m and see how the derivative on the plot changes
n = 2; % order of polynomial, reasonable numbers are 1-4
m = 5; % number of points in group to fit, m >= n+1
if m < n+1
    fprintf('you must increase m in the derivative polynomial fit! \n\n')
end
 
fprintf('see plot - hit return key to step \n')
 
% k = number of groups of m pts we can fit
% with dropping one pt and adding one pt for each group 
% Note this grouping is different than for quadrature
k = npts-m+1;
 
for j = 1:k
    i = j;
    % use Matlab std function polyfit to fit n-th
    % order polynomial to this group of pts
    % with polynomial coefficients returned in array a
    % with highest-order term first, lowest last
    a = polyfit(t(i:i+m-1),c(i:i+m-1),n);
        
    % differentiate polynomial analytically over interval
    % and evaluate at approx. mid-point time t(i+floor(m/2))
    sumd = 0; % mol/m3/s
    for jj = 1:n
        sumd = sumd + ... 
            a(jj)*(n-jj+1)*t(i+floor(m/2))^(n-jj);
    end
    
    % save derivative, time, and concentration values in arrays
    diffValue(j) = sumd; % mol/m3/s
    diffTime(j) = t(i+floor(m/2)); % s
    diffConc(j) = polyval(a,t(i+floor(m/2)));
    
    % OPTIONAL - SHOW PROCESS ON PLOT
    if i < 10 % just do this for first few groups
        % add the polynomial and the slope to the plot
        % get intercept b of the tangent line of the polynomial at the mid-pt
        b = polyval(a,t(i+floor(m/2))) - sumd*t(i+floor(m/2));
        % get start and end times
        tp(1) = t(i);
        tp(2) = t(i+m-1);
        % get y values of the tangent line at the start and end times
        yp = sumd*tp + b;
        % get conc values for the polynomial at each time
        cp = polyval(a,t(i:i+m-1));
        plot(t,c,'.')
        axis([0 1.1*max(t) 0 1.1*max(c)])
        title('HIT SPACE KEY TO STEP - differentiating data')
        ylabel('conc (mol/m3)')
        xlabel('time (s)')  
        hold on
        plot(t(i:i+m-1),cp,'b',tp,yp,'r')
        hold off
        pause
        % YOU NEED TO HIT ANY KEY TO CONTINUE STEPPING
        % SEE THE POLYNOMIALS FIT ON THE PLOT
    end
    % END OF OPTIONAL - SHOW PROCESS ON PLOT
end
 
% generate pts for a line at zero
zv = zeros(2);
zt = [t(1) t(npts)];
 
% make final plot
plot(t,c,'.',zt,zv,'g',diffTime,diffValue,'r')
title('pts = data, red line = derivative')
ylabel('conc, C (mol/m3) or dC/dt (mol/m3/s)')
xlabel('time (s)')  
 
pause
 
rate = -diffValue;
figure(1)
plot(log(diffConc), log(rate), 'bo')
 
% fit straight line through points in log-log plot
coef = polyfit(log(diffConc), log(rate),1);
% the coefficients of 1st order polynomial are the order and ln(k)
nfit = coef(1); % order
k325fit = exp(coef(2)); % k at 325
hold on
lnr = polyval(coef,log(diffConc)); % compute points for line through points
plot(log(diffConc),lnr,'k')
s = sprintf('ln(rate) vs. ln(conc), k = %4f, order = %4f',k325fit, nfit);
title(s)
hold off

% now show the model rate vs. the experimental rate estimates
figure(2)
ratefit = k325fit * diffConc .^ nfit;
plot(diffConc, rate, 'bo', diffConc, ratefit,'k')
title('rate vs. conc, circles are from dC/dt estimates, line = model rate')
axis([0 1.1*max(diffConc) 0 1.1*max(rate)])
 
E = 60e3; % J/mol
R = 8.3145; % J/K/mol
k300fit = k325fit * exp((E/R)*(1/325 - 1/300))
 
fprintf('DONE - see final plot \n')
 
% note the random error or scatter in these data
% which means we don't get the "true" values computed by the lab
%
% question for you: what does "true" mean in real world?
%
% order = 1.93 vs. "true" order = 2 in lab
%
% k325 = 8.84e-04 vs. "true" k325 = 6.36e-04 in lab
%
% k300fit =
%
%    1.3901e-04
%
% vs. k300 = 1.0e-04 "true" rate coefficient in lab for Ea = 60 kJ/mol