% Examples of integrating area under experimental data (quadrature) 
% using a variety of methods.
%
% SEE OUTPUT IN COMMAND WINDOW TO SEE AREA FROM EACH METHOD
% 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 BAR TO STEP - fitting polynomial to integrate under data')
ylabel('tracer conc (mol/m3)')
xlabel('time (s)')

[rows cols] = size(c); % size is std matlab function
npts = rows;

% dt = constant time interval between pts in this data set
% see time data above for first pt after t = 0
% other times are rounded
dt = 0.535; % s, time interval

%--------- RECTANGLES --------------

% just do rectangles...

sumq = 0; % mol-s/m3
for i = 1:npts-1
    sumq = sumq + dt*c(i);
end

% to get amt tracer injected (mol)
% need to multiply sumq (mol-s/m3) by 
% flow rate (m3/s)

amt = 1e-5*sumq; % mol = mol/s * mol-s/m3
fprintf('integrated moles by RECTANGLES = %8.6f \n',amt)

%--------- TRAPEZOIDS --------------

% now do trapezoids...

sumq = 0;
for i = 1:npts-1
    sumq = sumq + dt*(c(i) + c(i+1))/2;
end
    
amt = 1e-5*sumq; % mol = mol/s * mol-s/m3
fprintf('integrated moles by TRAPEZOIDS = %8.6f \n',amt)

%----- SIMPSON'S 1/3 RULE ----------

% NOTE Simpson's rules fit a polynomial to minimum number
% of points to get a unique fit: number pts = polynomial order + 1
% Since the fit is unique, Simpson determined the polynomial
% coefficients and integrated the polynomial for you already!

% now use Simpson's 1/3 rule
% (fit quadratic polynomial uniquely to groups of 3 equally spaced pts)  
% integral approx = sumq over i to npts-1 of
%  interval*(1/3)*(c(i) + 4*c(i+1) + c(i+2))
% where interval between pairs of pts here is 0.5 s

n = 2; % order of polynomial
m = n+1; % number of points in group to fit, m = n+1 for Simpson's

sumq = 0; % mol-s/m3
% k = number of groups of m pts we can fit
% floor is std matlab function that rounds down to nearest integer
k = floor((npts-1)/(m-1));
for j = 1:k
    i = 1 + (j-1)*(m-1);
    sumq = sumq + dt*(1/3)*(c(i) + 4*c(i+1) + c(i+2));
end

% we may have not fit some pts
p = mod(npts-1,k); % number of pts not fit
if p ~= 0
    fprintf('Simp 1/3 left off last %4i pts, catch with trapezoid \n',p)
    % last pt # of groups above was 1+k*(m-1)
    % very last pt # of all data is npts
    % catch these with trapezoid method
    for i = 1+k*(m-1):npts-1
        sumq = sumq + dt*(c(i) + c(i+1))/2;
    end
end

amt = 1e-5*sumq; % mol = mol/s * mol-s/m3
fprintf('integrated moles by SIMPSONS 1/3 = %8.6f \n',amt)

%----- SIMPSON'S 3/8 RULE ----------

% now use Simpson's 3/8 rule 
% (fit cubic polynomial uniquely to groups of 4 equally spaced pts) 
% integral approx = sumq over i to npts-1 of
%  interval*(3/8)*(c(i) + 3*c(i+1) + 3*c(i+2) + c(i+3))

n = 3; % order of polynomial
m = n+1; % number of points in group to fit, m = n+1 for Simpson's

sumq = 0; % mol-s/m3
% k = number of groups of m pts we can fit
% floor is std matlab function that rounds down to nearest integer
k = floor((npts-1)/(m-1));
for j = 1:k
    i = 1 + (j-1)*(m-1);
    sumq = sumq + ...
        dt*(3/8)*(c(i) + 3*c(i+1) + 3*c(i+2) + c(i+3));
end

% we may have not fit some pts
p = mod(npts-1,k); % number of pts not fit
if p ~= 0
    fprintf('Simp 3/8 left off last %4i pts, catch with trapezoid \n',p)
    % last pt # of groups above was 1+k*(m-1)
    % very last pt # of all data is npts
    % catch these with trapezoid method
    for i = 1+k*(m-1):npts-1
        sumq = sumq + dt*(c(i) + c(i+1))/2;
    end
end

amt = 1e-5*sumq; % mol = mol/s * mol-s/m3
fprintf('integrated moles by SIMPSONS 3/8 = %8.6f \n',amt)

% NOTE can get more efficient computation if pull
% dt*(3/8) outside of for-end repeat and multiple sumq
% by that amount after the for-end repeat but
% this is short and fast so there will be no noticeable
% speed difference here - similar for other repeats above

%----- GENERAL POLYNOMIAL FIT ----------

% We can also fit a polynomial to a greater number of
% points in each group to get some "smoothing" of the data
% Use a larger number of points (m) with more scatter in data
%
% ALSO, the following will also work with pts that are
% NOT equally spaced in time (Simpson's required equally spaced data)

n = 2; % order of polynomial, reasonable values are 1-3
m = 5; % number points in each group, m >= n+1
if m < n+1
    fprintf('you must increase m in the general polynomial fit! \n\n')
end

fprintf('see plot - hit return key to step \n')

sumq = 0; % mol-s/m3
% k = number of groups of m pts we can fit
% floor is std matlab function that rounds down to nearest integer
k = floor((npts-1)/(m-1));
for j = 1:k
    i = 1 + (j-1)*(m-1);
    % 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);
    
    % 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
    cp = polyval(a,t(i:i+m-1));
    hold on
    plot(t(i:i+m-1),cp,'r')
    hold off
    pause
    % YOU NEED TO HIT ANY KEY TO CONTINUE STEPPING
    % SEE THE POLYNOMIALS FIT ON THE PLOT
    % END OF OPTIONAL - SHOW PROCESS ON PLOT

    % integrate polynomial analytically over interval
    for jj = 1:n+1
        sumq = sumq + ... 
            (a(jj)/(n+2-jj))*(t(i+m-1)^(n+2-jj) - t(i)^(n+2-jj));
    end
end
    
% we may have not fit some pts
p = mod(npts-1,k);
if p ~= 0
    fprintf('Poly left off last %4i pts, catch with trapezoid \n',p)
    % last pt # of groups above was 1+k*(m-1)
    % very last pt # of all data is npts
    % catch these with trapezoid method
    for i = 1+k*(m-1):npts-1
        sumq = sumq + dt*(c(i) + c(i+1))/2;
    end
end

amt = 1e-5*sumq; % mol = mol/s * mol-s/m3
fprintf('integrated moles by POLYNOMIAL = %8.6f \n',amt)

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