function Cannon_d(d)
% d is the required travel distance of the cannonball
clc
fprintf('Required travel distance = %6.3f\n\n',d)
fprintf('Required initial velocity as a function of initial angle A (degrees) and headwind w\n\n')

fprintf(' A w = -20 w = -10 w = 0 w = 10 w = 20 \n')
fprintf('---------------------------------------------------------------------------')
for angle = [10 15 20 25 30 35 40 45 50 55 60 65 70]
fprintf('\n %2d ',angle)
for w = [-20 -10 0 10 20]
% For this initial velocity, how far does the cannonbal travel?
theta = pi*angle/180;
v0 = HowFast(theta,w,d);
fprintf(' %10.3f',v0)
end
end
fprintf('\n')

function vOpt = HowFast(theta,w,d)
vOpt = fzero(@(v0) F(v0,theta,d,w),[10 4000]);

function z = F(v0,theta,d,w)
% z is zero if the initial condition
u0 = [0 0 theta v0]';
opts = odeset('events', @g);
[t,uVals] = ode45(@(t,u) f(t,u,w), [0 Inf], u0, opts);
z = uVals(length(t),1)-d;

function up = f(t,u,w)
% t is a scalar (time in seconds)
% u is a column 4-vector with the property that
% u(1) = x(t), the x-coordinate of the cannonball at time t
% u(2) = y(t), the y-coordinate of the cannonball at time t
% u(3) = theta(t), the angle that the velocity vector makes with the horizontal
% u(4) = v(t), the speed of the projectile at time t
%
% w is the wind (+10 means tailwind of 10m/s, -10 means headwind of 10m/s)

% up is a column 4-vector that reports the rate of change of these 4
% functions at time t...
up = zeros(4,1);
g = 9.81; % acceleration of gravity
m = 15; % mass of the cannonball
c = 0.2; % the drag coefficient
rho = 1.29; % density of air
s = 0.25; % cross-sectional area of the cannonball

up(1) = u(4)*cos(u(3));
up(2) = u(4)*sin(u(3));
up(3) = -g/u(4)*cos(u(3));
Drag = (c*rho*s/2) * ( (up(1)+w).^2 + up(2).^2 );
up(4) = -Drag/m - g*sin(u(3));

function [gstop, isterminal, direction] = g(t,u,w)
% Event function to turn off ode45 as soon as the cannonball lands
gstop = u(2); % i.e., y(t) = 0
isterminal = 1;
direction = -1;

% reasonable bracketing interval guess with comments: 1pt
% putting together the right function for fzero: 3pt
% correct implementation: 1pt