- Published on

# Numerical Simulation of a Golf Shot - Part 3

- Authors
- Name
- Tinker Assist

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)
```