%% Shuttlecocks and the Akaike Information Criterion
%% Kurt Bryan and SIMIODE
% A Matlab script to help explore the project in Section 3.5.4.
%
% The Data: First, the data for the shuttlecock's fall, in time (seconds)/distance (meters) pairs:

times = [0, 0.347, 0.47, 0.519, 0.582, 0.650, 0.674, 0.717, 0.766, 0.823, 0.870, 1.031, 1.193, 1.354, 1.501, 1.726, 1.873];
dists = [0, 0.61, 1.00, 1.22, 1.52, 1.83, 2.00, 2.13, 2.44, 2.74, 3.00, 4.00, 5.00, 6.00, 7.00, 8.50, 9.50];

% A quick plot of (time, distance) pairs:

scatter(times,dists);

% The Model: We might posit a model of the form v'(t) = g (no air resistance) and consider g as an unknown, to be estimated. Then the governing ODE is (from equation (3.68) in the text)

syms v(t); %Declare v(t) as symbolic function
syms g;
ode = diff(v(t),t) == g %Define the ODE
vsol(t) = dsolve(ode,v(0)==0) %Incorporate initial condition

% Integrate to obtain the position in terms of g and t:

syms tau;
x(t,g) = int(vsol(tau),tau,0,t)

% Estimating Parameters: Form a sum of squares

syms SS(g)
SS(g) = sum((x(times,g)-dists).^2);

% To get a sense of where the minimum is, plot this expression as a function of k:

fplot(SS(g),[0 15])

% Solve SS'(g) = 0 to find the least-squares estimate for gravitational acceleration

eqn = diff(SS,g)==0;
gbest = vpasolve(eqn,g,[5 15])

% The residual is

SS(gbest)

% Use this value in X(t) to plot and compare to the data

fplot(x(t,gbest),[0 1.873],'-r')
hold on;
scatter(times,dists);
hold off;

%