% Class 21 - Numerical Integration
% Coded by NFR on 11.5.2019 
%
% Warm Up Example
% Another example of DoE experiment + optimization
% (should help with PSET 9)
%
D = readmatrix('Gluc_v2');
mdl = fitlm(D(:,1:2),D(:,3),'quadratic');
% Upon inspection we see that the x1:x2 cross term is high, so we will
% remove this, and refit the model
New_mdl = removeTerms(mdl,'x1:x2');
% REMINDER: only remove one term at a time
% We have a good model, we need to plot
% We need the model coeff/parameters
P = New_mdl.Coefficients.Estimate
% Define the model w/ these parameters
func = @(x,y) P(1) + P(2)*x + P(3)*y...
    + P(4)*x.^2 + P(5)*y.^2;
% Now we need to create matrix of x and y points
% Let's plot the full extent of our DoE variables
xvec = linspace(5,7,100);
yvec = linspace(1,2.5,100);
[xmat,ymat] = meshgrid(xvec,yvec);
zmat = func(xmat,ymat);
contourf(xvec,yvec,zmat,'ShowText','on')
xlabel('Growth Time')
ylabel('Exposure Time')
%
% Unit III - Numerical Integration
% Use the trapezoid method on the rocket data
% to determine the overal impulse. 
% What class of engine did we build?
D = readmatrix('Rocket_data');
nr = length(D);
Impulse = 0;
for i = 1:nr-1
    b = D(i+1,1); % Time point ahead
    a = D(i,1);  % Time point currently at
    Fa = D(i,2); % Thurst at a
    Fb = D(i+1,2); % Thrust at b
    Isection = (b-a)*(Fa+Fb)/2;
    Impulse = Impulse + Isection;
end
Impulse
% The other Newton Cotes Methods are similar
% They use more points to improve the accuracy
%
% Simpson's Rules (1/3 and 3/8 rule)
% 
% Book problem 19.4 (use simpson and others)
%
F1 = @(x) 1-x-4*x^3+2*x^5;
I1 = @(x) x - x^2/2-x^4+x^6/3; % Integrated form
a = -2;
b = 4;
% Part (a) analytical solution
Sol1 = I1(b)-I1(a)
fplot(F1,[a b])
% Part (b) single segement trapezoid
Sol2 = (b-a)/2*(F1(a)+F1(b))
% Part (c) Composite trapezoid - 2 segments
% NOTE: the coefficients don't change, 
% so we can calculate the composite in one line
Sol3 = (3)/2*(F1(a)+2*F1(1)+F1(b))
% for 4 segments (5 evaluation points!)
Sol4 = (1.5)/2*(F1(a)+2*F1(a+1.5)+2*F1(a+2*1.5)...
    +2*F1(a+3*1.5)+F1(b))
% Boole's rule (5 point)
% Table 19.2 - evaluating at 5 equdistant spaced points
% (same as above) except you will use
% different coefficients
%
% Matlab has a couple built in functions
% Practice these with Problem 19.8
%
t = [1 2 3.25 4.5 6 7 8 8.5 9 10];
v = [5 6 5.5 7 8.5 8 6 7 7 5];
% You could use a custom loop (see rocket example)
% or use the trapz function
Distance = trapz(t,v)
Avg_Velocity = Distance/9
% Part b - fit with a cubic equation
% then integrate that fit equation to find 
% the distance (show how we use polyint)
% Fit with cubic
p = polyfit(t,v,3);
% You can by hand determine the integrated form
% of this polynomial as we did before, or use polyint
pI = polyint(p);
Distance2 = polyval(pI,10)-polyval(pI,1)