function ballistics_distances
%BALLISTICS_DISTANCES: MATLAB function M-file that computes
%and plot the distances for all three cases, and also 
%computes the least squares error from the data in each 
%case.
%Author: Peter Howard
%
%Data
theta = 5:5:85;
dexp = [4.37 5.23 6.95 7.84 8.17 8.69 8.81 8.99 8.95 8.83 8.19 7.84 7.12 6.38 5.08 3.34 2.13];
%No air resistance
g = 9.81; H = .18; v = 10.26;
d1 = v*cosd(theta).*(v*sind(theta)+sqrt(v^2*sind(theta).^2 + 2*H*g))/g;
%
%Linear air resistance
b = .2792; v = 11.2343; 
%The implicit equation for distance d is
for k=1:17
    f = @(d) H + (g/b^2)*log(1-b*d/(v*cosd(theta(k)))) + (v*sind(theta(k))/b+g/b^2)*b*d/(v*cosd(theta(k)));
    guess = 10.26*cosd(theta(k))*(10.26*sind(theta(k))+sqrt(10.26^2*sind(theta(k))^2+2*H*g))/g;
    %I.e., why not use the distance associated with no air resistance as an 
    %initial approximation. 
    d2(k) = fzero(f,guess);
end
%
%Nonlinear air resistance
for k=1:17
    d3(k) = dist(theta(k));
end
%Plots and errors
plot(theta,dexp,'o',theta,d1,theta,d2,'--',theta,d3,':')
title('Plot of Horizontal Distances, model and data','FontSize',15)
xlabel('Angle of Inclination')
ylabel('Horizontal Distance Traveled')
legend('Experimental data','No air resistance','Linear air resistance','Nonlinear air resistance')
error1 = norm(dexp-d1)
error2 = norm(dexp-d2)
error3 = norm(dexp-d3)
%
function d = dist(theta)
%DIST: MATLAB function M-file that takes an angle theta as input
%and returns the horizontal distance a dart with initial speed
%v0 and air resistance coefficient b
v0 = 12.2540; b=.0645; g = 9.81;
rhs = @(t,y) [y(2);-b*sqrt(y(2)^2+y(4)^2)*y(2);y(4);-g-b*sqrt(y(2)^2+y(4)^2)*y(4)];
options=odeset('Events',@subevents);
[t,y,te,ye,ie]=ode45(rhs,[0 20],[0;v0*cosd(theta);.18;v0*sind(theta)],options);
d = ye(end,1);
%
function [lookfor stop direction]=subevents(t,y)
lookfor = y(3);
stop = 1;
direction = -1;