#Shuttlecocks and the Akaike Information Criterion

#A worksheet 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:

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

N = len(shuttledata);

#A plot
plt1 = scatter_plot(shuttledata,facecolor='red');
show(plt1)

#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)
var('g t');
v = function('v')(t)
de = diff(v,t) == g
vsol(t) = desolve(de,v,[0,0],t); #Solve with v(0) = 0.

#Integrate to obtain the position in terms of t and k, taking x(0) = 0.
var('tau');
x(t,g) = integral(vsol(tau),tau,0,t);

#Estimating Parameters: Form a sum of squares
SS = function('SS')(g);
SS(g) = add((x(shuttledata[i][0])-shuttledata[i][1])^2 for i in range(N));

#A quick plot of SS(g)
plot(SS(g), (g,0,15)).show()

#We can solve by setting dSS/dg = 0 and solving for g.
criteqn = diff(SS,g) == 0;
gbest = find_root(criteqn,5,15) #Look between 5 and 15.

#The residual is
SS(gbest)

#Plot the function x(t) with this best choice for g.
plt2 = plot(x(t,gbest), t, 0, 1.873);
pp = plt1 + plt2;
show(plt1+plt2)

#An alternative is to use Sage's "model fitting" ability, to adjust
#k in the function x(t,g) to best fit the data.
model(t) = x(t,g); #Specify the model to be fit, with "t" as the independent variable
sol = find_fit(shuttledata,model,parameters=[g]) #Fit the model by adjusting k and b
bestx(t) = model(g=sol[0].rhs()) #Define the best fit function of the form in "model"
bestx(t)