function golf_ball_solutions

% Model of the flight of a golf ball for various lift and drag
% considerations. Uses 'ode45' as the numerical integrator. 
%
% Joe Schoneman - 20 Jan 13 

clc;clear;close all % Clear the command window, workspace, and 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]
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 velocities                         
time_interval      = [0 t_span];   % Interval of solution

% Solution 1: No lift, no drag
Cd = 0; Cl = 0; % Coefficients of lift and drag
Kd = Cd*rho*A/2; Kl = Cl*rho*A/2; % Lumped lift and drag coefficients
[t1 w1] = ode45(@golf_ode, time_interval, initial_conditions);

x1 = w1(:,1); % Extract x and y values from the w vector
y1 = w1(:,2);
x1 = x1(y1 >= 0); % Logical indexing: only interested in positive y-values
y1 = y1(y1 >= 0);

% Solution 2: No lift, drag
Cd = .23; Cl = 0; % Coefficients of lift and drag
Kd = Cd*rho*A/2; Kl = Cl*rho*A/2; % Lumped lift and drag coefficients
[t2 w2] = ode45(@golf_ode, time_interval, initial_conditions);

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

% Solution 3: Lift and drag
Cd = .23; Cl = .12; % Coefficients of lift and drag
Kd = Cd*rho*A/2; Kl = Cl*rho*A/2; % Lumped lift and drag coefficients
[t3 w3] = ode45(@golf_ode, time_interval, initial_conditions);

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

% Now, extract x and y values from the w vector, as we assigned below.


% Analytical Solution:

x_a = 0.1:500;
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 all solutions:

% - Analytical solution as a blue line;
% - Numerical with no drag/ift as green circles
% - Numerical with drag/no lift as red circles
% - Numerical with drag and lift as blue circles

% Plot Setup (Use Plot Editor to auto-generate code like this)

figure1 = figure;
axes1 = axes('Parent',figure1);
box(axes1,'on');
hold(axes1,'all');

plot(x_a,y_a,'Parent',axes1,'DisplayName','Analytical');
plot(x1,y1,'Parent',axes1,'Marker','o','LineStyle','none','Color',[1 0 1],...
    'DisplayName','Numerical, No Lift/Drag');
plot(x2,y2,'Parent',axes1,'Marker','o','LineStyle','none','Color',[1 0 0],...
    'DisplayName','Numerical, Drag Only');
plot(x3,y3,'Parent',axes1,'Marker','o','LineStyle','none','Color',[0 0 1],...
    'DisplayName','Numerical, Lift and Drag');
xlabel('x [m]');
ylabel('y [m]');
title('Analytical vs Various Numerical Trajectories');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.271620133399293 0.137748896807607 0.323798093911217 0.161441616482855]);


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

function eom = golf_ode(t,w)

global m Kl Kd g

% Assign values from input vector 'w'(row vector)

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]; % Output Vector 'eom' (Column vector)