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