If you have yet read parts 1 and 2, I recommend checking them out to see how we got here

Golf Shot Simulation - Part 1 - Golf Ball Aerodynamics

Golf Shot Simulation - Part 2 - Numerical Integration Methods

Define Simulation Function

We'll start by defining a function "golfshot_sim", where we will run the full simulation

def golfshot_sim(initial_velocity, launch_angle, rotation_rate):
First, we need to define the constants we will be referencing

g = -9.81 # acceleration due to gravity (m/s/s)
r = 0.0215 # radius of a golf ball (m)
m = 0.0459 # mass of a golf ball (kg)
D = r*2 # diameter of a golf ball
A = math.pi*r**2 # cross sectional area of a golf ball
kin_visc = 1.48*10**-5 # kinematic viscosity of air at 15*C ( m^2/s )
rho = 1.225 # density of air at sea level, 15*C, no humidity ( kg/(m^3) )

Now we need to define initial state vectors for position, velocity, and acceleration

# initial position is our origin - all zeros
pos0 = [0, 0, 0]

# initial velocity vector is a function of our launch angle and launch velocity magnitude
vel0 = [initial_velocity*math.cos(math.radians(launch_angle)),
        initial_velocity*math.sin(math.radians(launch_angle)),
        0]

# only initial acceleration is gravity
accel0 = [0, -g, 0]

# combine acceleration, velocity, and position into numpy state array
state = np.array([[pos0, vel0, accel0]])

Define numerical integration loop parameters such as the simulation length in time, and total steps. From those values, we can compute the time of each step size and the full time vector.

# parameters for executing the simulation
sim_time = 15 # 100s should be more than plenty of time
steps = 30000 # number of steps in the simulation
delta_t = sim_time/steps # change in time between each sim step
time = np.linspace(0, sim_time, num=steps) # time vector for sim
Now we must define our simulation loop, where each state will be computed.

# start the numerical simulation
for curr_time in time:
    # compute magnitude of our velocity vector
    V = np.linalg.norm(state[-1,1,:])

    # get reynolds number for the current ball velocity
    Re = V*D/kin_visc

    # Compute drag coefficient based on reynolds numbers. Diff equations for high and low speeds
    if Re < 10**5:
        c_D = 1.29*10**-10*Re**2 - 2.59*10**-5*Re + 1.50   
    else:
        c_D = 1.91*10**-11*Re**2 - 5.40*10**-6*Re + 0.56 

    # compute acceleration due to drag on the golf ball
    a_D = V**2*c_D*rho*A/(2*m)
    
    # compute non-dimensional spin factor S
    S = rotation_rate*r/V

    # compute the lift coefficient. Function of spin factor
    c_L = -3.25*S**2 + 1.99*S

    # compute acceleration due to lift
    a_L = V**2*c_L*rho*A/2/m 

    # compute current ball path angle from velocity vector
    angle = math.atan(state[-1,1,1]/state[-1,1,0]) 

    # compute acceleration in each axis
    # drag acceleration acts opposite the ball path
    # lift acceleration acts orthogonal to ball path
    a_x = -a_D*math.cos(angle) - a_L*math.sin(angle)
    a_y = -a_D*math.sin(angle) + a_L*math.cos(angle) - abs(g)
    a_z = 0

    # compute velocity in each axis
    v_x = state[-1,1,0] + state[-1,2,0]*delta_t
    v_y = state[-1,1,1] + state[-1,2,1]*delta_t
    v_z = 0

    # compute position
    x = state[-1,0,0] + state[-1,1,0]*delta_t
    y = state[-1,0,1] + state[-1,1,1]*delta_t
    z = 0
    
    # consolidate state array and append to previous states
    state = np.append(state,[[[x, y, z], [v_x, v_y, v_z], [a_x, a_y, a_z]]],0)

    # if the height of the ball drops below 0 on the y axis, we know the ball has hit the ground
    if y < 0:
        break

Now that we have our full state history as a numpy array, we will convert this to a pandas timeseries and output that data structure

# pull out x and y states for 2D sim
pos_x = (state[:,0,0]).tolist()
pos_y = (state[:,0,1]).tolist()

# shorten our time vector to include only the values between 0 and time when ball crossed y-axis
time = time[0:len(pos_x)]

# Create pandas timeseries dataframe
np_2d_timeseries = np.array([time, pos_x, pos_y])
pd_2d_timeseries = pd.DataFrame(np_2d_timeseries.T, columns=['time','x','y'])

# return pandas timeseries dataframe
return pd_2d_timeseries

Running the Simulation

Now that we have defined our function, let's pass in some initial values and run it

state_history = golfshot_sim(launch_velocity=59,launch_angle=12.6,rotation_rate=342)
print(state_history)
And below is an example of what is returned in the terminal

time           x         y
0      0.000000    0.000000  0.000000
1      0.000500    0.028790  0.006435
2      0.001000    0.057579  0.012873
3      0.001500    0.086365  0.019311
4      0.002000    0.115146  0.025748
...         ...         ...       ...
14566  7.283243  230.060218  0.031105
14567  7.283743  230.070649  0.022415
14568  7.284243  230.081079  0.013724
14569  7.284743  230.091510  0.005032
14570  7.285243  230.101940 -0.003661

[14571 rows x 3 columns]

golfshot Python Package

Instead of writing all of this code out on your own, you can install the golfshot python package.

> pip install golfshot

from golfshot import GolfShot

sim = GolfShot(launch_velocity=75,launch_angle=10,rotation_rate=300)
state_history = sim.run_sim()

print(state_history)

Full Code

import math
import numpy as np
import pandas as pd

def golfshot_sim(launch_velocity, launch_angle, rotation_rate):
    g = -9.81 # acceleration due to gravity (m/s/s)
    r = 0.0215 # radius of a golf ball (m)
    m = 0.0459 # mass of a golf ball (kg)
    D = r*2 # diameter of a golf ball
    A = math.pi*r**2 # cross sectional area of a golf ball
    kin_visc = 1.48*10**-5 # kinematic viscosity of air at 15*C ( m^2/s )
    rho = 1.225 # density of air at sea level, 15*C, no humidity ( kg/(m^3) )
    
    # initialize state - position, velocity, and acceleration
    # simulation is 2D. Thus, all z componenents are zero

    # initial position is our origin - all zeros
    pos0 = [0, 0, 0]

    # initial velocity vector is a function of our launch angle and launch velocity magnitude
    vel0 = [launch_velocity*math.cos(math.radians(launch_angle)),
            launch_velocity*math.sin(math.radians(launch_angle)),
            0]

    # only initial acceleration is gravity
    accel0 = [0, -g, 0]

    # combine acceleration, velocity, and position into state array
    state = np.array([[pos0, vel0, accel0]])

    # parameters for executing the simulation
    sim_time = 15 # 100s should be more than plenty of time
    steps = 30000 # number of steps in the simulation
    delta_t = sim_time/steps # change in time between each sim step
    time = np.linspace(0, sim_time, num=steps) # time vector for sim

    # start the numerical simulation
    for curr_time in time:
        # compute magnitude of our velocity vector
        V = np.linalg.norm(state[-1,1,:])

        # get reynolds number for the current ball velocity
        Re = V*D/kin_visc

        # Compute drag coefficient based on reynolds numbers. Diff equations for high and low speeds
        if Re < 10**5:
            c_D = 1.29*10**-10*Re**2 - 2.59*10**-5*Re + 1.50   
        else:
            c_D = 1.91*10**-11*Re**2 - 5.40*10**-6*Re + 0.56 

        # compute acceleration due to drag on the golf ball
        a_D = V**2*c_D*rho*A/(2*m)
        
        # compute non-dimensional spin factor S
        S = rotation_rate*r/V

        # compute the lift coefficient. Function of spin factor
        c_L = -3.25*S**2 + 1.99*S

        # compute acceleration due to lift
        a_L = V**2*c_L*rho*A/2/m 

        # compute current ball path angle from velocity vector
        angle = math.atan(state[-1,1,1]/state[-1,1,0]) 

        # compute acceleration in each axis
        # drag acceleration acts opposite the ball path
        # lift acceleration acts orthogonal to ball path
        a_x = -a_D*math.cos(angle) - a_L*math.sin(angle)
        a_y = -a_D*math.sin(angle) + a_L*math.cos(angle) - abs(g)
        a_z = 0

        # compute velocity in each axis
        v_x = state[-1,1,0] + state[-1,2,0]*delta_t
        v_y = state[-1,1,1] + state[-1,2,1]*delta_t
        v_z = 0

        # compute position
        x = state[-1,0,0] + state[-1,1,0]*delta_t
        y = state[-1,0,1] + state[-1,1,1]*delta_t
        z = 0
        
        # consolidate state array and append to previous states
        state = np.append(state,[[[x, y, z], [v_x, v_y, v_z], [a_x, a_y, a_z]]],0)

        # if the height of the ball drops below 0 on the y axis, we know the ball has hit the ground
        if y < 0:
            break

    # pull out x and y states for 2D sim
    pos_x = (state[:,0,0]).tolist()
    pos_y = (state[:,0,1]).tolist()

    # shorten our time vector to include only the values between 0 and time when ball crossed y-axis
    time = time[0:len(pos_x)]

    # Create pandas timeseries dataframe
    np_2d_timeseries = np.array([time, pos_x, pos_y])
    pd_2d_timeseries = pd.DataFrame(np_2d_timeseries.T, columns=['time','x','y'])

    # return pandas timeseries dataframe
    return pd_2d_timeseries
    
state_history = golfshot_sim(launch_velocity=59,launch_angle=12.6,rotation_rate=342)
print(state_history)