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)