% Class 22 - Integration Part II + Numerical Derivative
%
% First method (combines multiple estimates)
% to take integral from function
% Use Romberg to combine low quality trap estimates
% to form a higher accuracy estimate
%
% Problem 20.23
m = 80; %kg
Cd = 0.2; %kg/m
g = 9.8; %m/s^2
vfunc = @(t) sqrt(g*m/Cd)*tanh(sqrt(g*Cd/m)*t);
Distance = romberg(vfunc,0,8,1)
%
% Second Method - strategically picks points
% at which to evaluate your function
% Gauss Quadrature
% We can use table 20.1 after change of variables
% Problem 20.3
f1 = @(xd) (1.5+1.5*xd)*exp(2*(1.5+1.5*xd))*1.5;
co = 1;
c1 = 1;
xo = -1/sqrt(3);
x1 = 1/sqrt(3);
Int = co*f1(xo)+c1*f1(x1)
% Check answer with Romberg
% use the regular function
f2 = @(x) x.*exp(2*x);
Int2 = romberg(f2,0,3,0.5)
% You can also use built in functions
% Adaptive quadrature takes into account
% the form of the function
Int3 = quadgk(f2,0,3);
% More common built in functuion
% (Dr. Reuel's go-to)
Int4 = integral(f2,0,3)
% Also note there are built in functions
% for higher dimension - integral2 and integral3
%
% Numerical Derivative
% 1) Romberg to combine two poor estimates
% 2) Noisy (real) data should have a model
%   fit first before finding derivative
% 3) Built in matlab functions
% Return to the rocket problem
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];
% task is to plot acceleration vs. time
diff_t = diff(t)
diff_v = diff (v)
Accel = diff_v./diff_t;
% Note, we have one less one timepoint
% we should plot these values at the midpoints
n = length(t);
tm = (t(1:n-1)+t(2:n))./2;
plot(tm,Accel,'o')
xlabel('Time (s)')
ylabel('Accel (m/s^2)')
% Another built in function = gradient
% Demo with equat spacings
Tvec = [24 15 10 9 6 3 2];
% These temp values collected every five minutes
dTemp_dt = gradient(Tvec,5)
% Also note, this can be used for unequal spaced points
% except h would be a vector of distances
% Try with rocket problem
% Should be able to use h vector, will post update