% 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

% The experiment here involved injection of a pulse of
% 1e-3 mol of tracer into fluid flowing at 1e-5 m3/s
% through a packed bed.
% This is a Residence Time Distribution (RTD) experiment
% that is used to investigate flow patterns in process units.
% These data are time and tracer concentration measurements
% at the outlet of the bed.
% Data pasted in from experiment in SimzLab > ReactorLab >
% Division 4 Lab 2 (www.SimzLab.com)

% paste header data from experiment and comment out:

% D4L2, Data Set 1
% 
% input_1 remained constant at pulse 
% dp_1 remained constant at    2.00E-3  m
% wt_1 remained constant at      0.250  kg
% flow_1 remained constant at    1.00E-5  m3/s
% dBed_1 remained constant at      0.200  m
% len_1 remained constant at    1.14E-2  m
% t_1 has units of s
% Cout_1 has units of mol/m3
% 
% t_1	Cout_1

% paste data from experiment into array d
% alternatively, could also read data from a text file
% but do this here to keep everything in one file

d = [
-2.14	5.61E-2
-1.60	-0.157
-1.07	0.298
-0.535	6.79E-2
0	0.366
0.535	-0.139
1.07	0.242
1.60	0.352
2.14	1.65
2.67	2.10
3.21	3.19
3.75	4.02
4.28	5.17
4.82	5.82
5.35	6.54
5.88	7.05
6.42	7.78
6.96	8.15
7.49	8.05
8.03	8.42
8.56	8.31
9.10	7.83
9.63	7.74
10.2	7.37
10.7	7.41
11.2	7.03
11.8	6.54
12.3	6.26
12.8	5.54
13.4	5.24
13.9	4.94
14.4	4.67
15.0	4.33
15.5	3.49
16.1	3.49
16.6	3.09
17.1	3.01
17.6	2.42
18.2	2.28
18.7	1.73
19.3	1.82
19.8	1.29
20.3	1.41
20.9	0.915
21.4	0.812
21.9	1.11
22.5	0.860
23.0	0.779
23.5	0.826
24.1	0.606
24.6	0.396
25.1	0.268
25.7	0.601
26.2	0.356
26.8	0.200
27.3	0.174
27.8	0.341
28.4	-9.55E-2
28.9	-0.244
29.4	9.45E-2
30.0	0.136
30.5	6.72E-2
31.0	-4.32E-2
31.6	0.211
32.1	-8.75E-2
];

% extract time and concentration columns into separate arrays
t = d(:,1); % s, time
c = d(:,2); % mol/m3, tracer concentration

% plot
plot(t,c,'.')
title('HIT SPACE KEY TO STEP - differentiating data')
ylabel('tracer 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
    % >>> FREEMAT <<< SEEMS TO HAVE PROBLEMS WITH THIS
    % SECTION SO BEST TO COMMENT OUT THIS SECTION
    % IN FREEMAT UNDER TOOLS MENU
    if i < 10 % just do this for first 10 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,'.')
        title('HIT SPACE KEY TO STEP - differentiating data')
        ylabel('tracer 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,'.',diffTime,diffValue,'r',zt,zv,'g')
title('RTD experiment: pts = data, red line = derivative')
ylabel('tracer conc, C (mol/m3) or dC/dt (mol/m3/s)')
xlabel('time (s)')  

% FOR SOME APPLICATIONS such as the differential method
% of reaction kinetic analysis, you may need to find and delete
% (or find and keep) derivative values that are positive (or negative)
% before you take the logarithm of the values. This can be done using
% Matlab's find command, e.g.,
%     y = find(diffValue < 0); % y is array of indexes of diffValue where values are < 0
%     diffVx = diffValue(y); % diffVx is array with only the selected values of diffValue

fprintf('DONE - see final plot \n')