function elev_offset=softball_function(U,theta0)
% Perform ode45 solve of the equations of motion given in projectile.m.
%
% U - launch speed
% theta0 - launch angle (relative to the horizontal)
% elev_offset - difference between max(elevation) based on the ODE solve 
%               and ypeak, which is specified in the problem statement
 
global ypeak

% Specify maximum integration time.
tmax = 10.0;
 
% Solve ODEs using ode45 taking care to terminate the integration once v 
% becomes negative. The matrix y consists of a column vector for x 
% (horizontal position), a column vector for u (horizontal velocity), a
% column vector for y (elevation) and a column vector for v (vertical 
% velocity).
options = odeset('Events',@vneg,'MaxStep',1e-3);
[~,y] = ode45(@projectile,[0.0 tmax],[0.0 U*cos(theta0) ...
              0.0 U*sin(theta0)],options);

elev_offset = max(y(:,3))-ypeak;

%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy=projectile(~,y)
% The equations of motion for a projectile. 
% dy(1)=x
% dy(2)=u
% dy(3)=y
% dy(4)=v

global Cd Uwind g

dy = zeros(4,1); % Initialize.

dy(1) = y(2);
dy(2) = -Cd*(y(2)*sqrt(y(2)^2.0+y(4)^2.0)+Uwind^2.0);
dy(3) = y(4);
dy(4) = -Cd*y(4)*sqrt(y(2)^2.0+y(4)^2.0)-g;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [value,isterminal,direction]=vneg(~,y)

value = y(4);      % Focus on v, not u, x or y.
isterminal = 1;    % Stop ODE solve when v<0.
direction = -1;    % Only detect a decreasing crossing.