function WarmupQuiz
%Coded by NEK
close all

%Get data from file
%Data is in meters with 1s time intervals
Altitudes = xlsread('AltitudeData.xlsx');

%Time interval -- this is more important when it isn't 1
Interval = 1;

%Differentiate with respect to time
Acceleration = diff(Altitudes,2)/Interval;

%Plot the data
plot(Acceleration/9.8);
xlabel('Time, s')
ylabel('Acceleration, xG')

end

function ReuelRocketStuff
%Coded by N.K.
close all

%Get Data
DATA = xlsread('Rocket_data.xls');
%Time, s
T = DATA(:,1);

%Force, N
F = DATA(:,2);

%Plot Force vs Time
figure(1)
subplot(1,2,1)
hold on
plot(T, F);
xlabel('Time, s')
ylabel('Force, N')

%Time intervals are variable
subplot(1,2,2)
hold on
plot(T,cumtrapz(T,F));
xlabel('Time, s')
ylabel('Impulse, N')

%Total Impulse Trapezoidal
Impulse = trapz(T,F)


%If we wanted to use Newton-Cotes with one step, order 5
%This method would be best avoided, especially at high powers, where error
%increases exponentially.  If the order of the polyfit function is raised
%to 6, you can see error explode
P = polyfit(T,F,5);
Xmodel = linspace(0,T(end),100);
Ymodel = polyval(P,Xmodel);
subplot(1,2,1)
plot(Xmodel,Ymodel, '--')
subplot(1,2,2)
plot(Xmodel, cumtrapz(Xmodel,Ymodel), '--')
legend('Trapezoidal Integration', 'Fitted Polynomial')

%We can get an exact integral of the polynomial function.  If we were to
%add more steps, we would just make this a piecewise function and integrate
%each piece
Q = polyint(P);
ImpulseNC = polyval(Q,T(end))-polyval(Q,T(1))
end

function IntegratePolynomialwithSimpsonThing
%f(x) = 0.2 + 25x -200x^2 + 675x^3 - 900x^4 + 400x^5
%integrate from 0 to 0.8

%Polynomial
P = transpose([400; -900; 675; -200; 25; 0.2]);

%number of steps
n = 4;

%Bounds
a=0;
b = 0.8;

%X-values
X = linspace(a,b,n+1);

%Yvalues
Y = polyval(P,X);

%This function integrates with the 1/3 rule.  It will execute with an odd
%number of steps, but its output will be erroneous
Summation = Y(1);
%Sum the weighted y-values within a loop
for i = 2:n
    %a modulus function (mod(a,b) = remainder after dividing a by b) can be
    %used with an iterator to do slightly different operations every b'th
    %iteration
    %The NOT function (~) can be used with this to make it more robust
    if mod(i,2)
        Summation = Summation + 2*Y(i);
    else
        Summation = Summation + 4*Y(i);
    end
end
Summation = Summation + Y(end);


I = Summation*(b-a)/(3*n)

%Calculate 4th derivative
%We're cheating here by calculating the actual derivative.  We could 
%calculate this numerically, but we only have 5 points total and would be
%left with only one point in the fourth derivative
Y_4 = polyder(polyder(polyder(polyder(P))));
Y_4_avg = (polyval(polyint(Y_4),b)-polyval(polyint(Y_4),a))/(b-a);

%Calculate the Error
EstimatedError = -(b-a)^5/(180*n^4)*Y_4_avg

%Exact error
Q = polyint(P);
ActualError = I - diff(polyval(Q, [a b]))
end