function golf_ball_example

% Model of the flight of a golf ball for comparison of analytical model with
% a basic lift/drag model. Uses 'ode45' as the numerical integrator. 

clc;clear;close all % Clear the command window, workspace, and close any figures

global m Kl Kd g % Variables that must be present, but are not passed, to 
                 % the ODE function

% Some starting parameters of interest

m = 45.93e-3; % Ball mass [kg]
d = 42.67e-3; % Ball diameter [m]
V_0 = 165; % Starting velocity [Mi/hr]
theta_0 = 12; % Launch angle [deg]
rho = 1.225; % Atmospheric density [kg/m^3]; %
g = 9.81; % Gravitational acceleration [m/s]

% Derived Parameters

V_0 = V_0 * 0.44704; % Starting velocity [m/s]
theta_0 = theta_0 * pi/180; % Lauch angle [rad]
A = d^2/4; % Cross-sectional ball area [m^2]

% Initial conditions 
Vx = V_0*cos(theta_0); % Velocity x-component
Vy = V_0*sin(theta_0); % Velocity y-component
t_span = 4.5; % Estimate time in flight
initial_conditions = [0 0 Vx Vy];  % Initial conditions for integrator                     
time_interval      = [0 t_span];   % Time inteval for integrator


% Numerical Solution

Cd = .23; Cl = .12; % Coefficients of lift and drag
Kd = Cd*rho*A/2; Kl = Cl*rho*A/2; % Lumped lift and drag coefficients
[t w] = ode45(@golf_ode, time_interval, initial_conditions); % Call solution

x = w(:,1); % Extract x and y values from the w vector
y = w(:,2);
y = y(y > 0); % Only interested in positive y-values
x = x(y > 0); 

% Analytical Solution:

x_a = 0.1:250;
y_a = -g/2*(sec(theta_0)/V_0)^2 * x_a.^2 + tan(theta_0)*x_a; % Parabolic solution
x_a = x_a(y_a > 0);
y_a = y_a(y_a > 0); 

% Plot Results:

plot(x,y,'or',x_a,y_a); % 'or' means 'circles, red'
xlabel('x [m]');
ylabel('y [m]');
title('Analytical vs Numerical Trajectories');

%%%%%%%% END MAIN FUNCTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function eom = golf_ode(t,w)

global m Kl Kd g

% Assign values from input vector {w} (What We Want)

x = w(1);
y = w(2);
xdot = w(3); % dx/dt
ydot = w(4); % dy/dt

e1 = xdot;
e2 = ydot;
e3 = 1/m * sqrt(xdot^2 + ydot^2) * (-Kd*xdot - Kl*ydot); % EOM for x
e4 = 1/m * sqrt(xdot^2 + ydot^2) * (Kl*xdot - Kd*ydot) - g; % EOM for y

eom = [e1; e2; e3; e4]; % {eom} vector