function ballistics_trajectories
%BALLISTICS_TRAJECTORIES: MATLAB function M-file for plotting
%trajectories at 35 degrees associated with the ballistics 
%project for the cases of no air resistance, linear air resistance, 
%and nonlinear air resistance. 
%Author: Peter Howard
%
%No air resistance
g = 9.81; H = .18; theta = 35; v = 10.26;
tground1 = (v*sind(35)+sqrt(v^2*sind(35)^2+2*H*g))/g;
t1 = linspace(0,tground1,50);
x1 = v*cosd(theta)*t1;
y1 = -g*t1.^2/2+v*sind(theta)*t1+H;
%
%Linear air resistance
b = .2792; v = 11.2343;
f = @(t) H - g*t/b + (v*sind(theta)/b+g/b^2)*(1-exp(-b*t));
tground2 = fzero(f,tground1);
t2 = linspace(0,tground2,50);
x2 = (v/b)*cosd(theta)*(1-exp(-b*t2));
y2 = H - g*t2/b + (v*sind(theta)/b+g/b^2)*(1-exp(-b*t2));
%
%Nonlinear air resistance
v = 12.2540; b=.0645; 
rhs = @(t3,y3) [y3(2);-b*sqrt(y3(2)^2+y3(4)^2)*y3(2);y3(4);-g-b*sqrt(y3(2)^2+y3(4)^2)*y3(4)];
options=odeset('Events',@subevents);
[t3,y3,te,ye,ie]=ode45(rhs,linspace(0,3,1000),[0;v*cosd(theta);.18;v*sind(theta)],options);
%
plot(x1,y1,x2,y2,y3(:,1),y3(:,3))
axis equal;
axis([0 12 0 3]);
title('Trajectories for \theta = 35^o');
legend('No air resistance', 'Linear air resistance', 'Nonlinear air resistance');
xlabel('Horizontal distance')
ylabel('Height')
%
function [lookfor stop direction]=subevents(t,y3)
lookfor = y3(3);
stop = 1;
direction = -1;