"""golf.py: a module for solving the golf ball projectile motion problem
Copyright (C) 2009-2010 Thomas J. Duck
All rights reserved.
Thomas J. Duck
Department of Physics and Atmospheric Science,
Dalhousie University, Halifax, Nova Scotia, Canada, B3H 3J5
NOTICE
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License version 2 as
published by the Free Software Foundation.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston,
MA 02111-1307 USA
"""
import numpy
import numpy.random; numpy.random.seed()
import scipy.optimize
from numpy import pi
h = 0.01 # Time step (s)
g = 9.8 # Gravity (m/s^2)
dens = 1.2 # Air density (kg/m^3)
m = 45.93e-3 # Mass of golf ball (kg)
d = 42.67e-3 # Ball diameter (m)
A = pi*(d/2.)**2 # Cross-sectional area of golf ball (m^2)
def set_params(h_=h, g_=g, dens_=dens, m_=m, d_=d, A_=A):
global h, g, dens, m, d, A
h,g,dens,m,d,A = h_,g_,dens_,m_,d_,A_
def get_spin_ratio(omega, V):
"""Returns the spin ratio.
A ValueError is raised if the spin ratio is outside the range allowed
by the lift and drag parameterizations.
"""
R = omega*(d/2.)/V
if R<0.05 or R>3:
raise ValueError, 'Spin ratio outside of allowed range: %f' % R
else:
return R
def get_Cd(omega, V):
"""Returns the drag coefficient Cd
The parameterization from the USGA Second Report on the Study of
Spin Generation (2007, Appendix C) is used.
"""
R = get_spin_ratio(omega, V)
return 0.1403 - 0.3406*R*numpy.log(R) + 0.3747*R**1.5
def get_Cl(omega, V):
"""Returns the lift coefficient Cl.
The parameterization from the USGA Second Report on the Study of
Spin Generation (2007, Appendix C) is used.
"""
R = get_spin_ratio(omega, V)
return 0.3996 + 0.1583 * numpy.log(R) + 0.03790 * R**(-0.5)
def get_U_default(x,z,t):
return 0.
def get_zground_default(x):
return get_zground_default.f(x)
get_zground_default.f = lambda x: 0.
# Perform the finite difference calculation
def get_trajectory(V0, theta0, rpm0, x0=0, z0=0,
get_zground=get_zground_default,
get_U=get_U_default,
N=100, t0=0.):
"""Returns the trajectory x,z coordinates for a golf ball.
V0 - the launch speed (m/s)
theta0 - the launch angle (degrees)
rpm0 - the rotation rate of the ball at launch (revolutions per minute)
x0 - initial horizontal distance (m)
z0 - initial altitude (m)
get_zground(x) - a function that returns the ground altitude (m) for
position x
get_U(x,z,t) - a function that returns the eastward wind speed (m/s) for
position x,z and time t
N - the maximum number of points to return in the trajectory
t0 - the time (s) at the beginning of the trajectory
The calculation is continued until the projectile is below ground.
Although some points may not be saved along the trajectory, the last
two points are guaranteed to be neighbours so that landing and speed
calculations may be performed.
"""
raise RuntimeError, "Don't use this; use cgolf.get_trajectory() instead"
theta0 = numpy.radians(theta0) # Launch angle (radians)
u0 = V0 * numpy.cos(theta0) # Initial horizontal ground speed (m/s)
w0 = V0 * numpy.sin(theta0) # Initial vertical ground speed (m/s)
# Initialize previous and current values for x, z and V (net ball air speed)
xprev, xcur = x0, u0*h + x0
zprev, zcur = z0, w0*h + z0
Vprev, Vcur = numpy.sqrt((u0 - get_U(x0, z0, t0))**2 + w0**2), None
# Create arrays to store the trajectory, and save the first points
x = numpy.zeros(N) # Horizontal distances (m)
z = numpy.zeros(N) # Altitudes (m)
filled = numpy.zeros(N,dtype=numpy.bool) # Flags points filled in
x[0], z[0], filled[0] = xprev, zprev, True
x[1], z[1], filled[1] = xcur, zcur, True
# Determine the number of steps between save points. Assume flight
# times are up to 15 seconds long.
step = int(numpy.ceil(round(15./h)/N))
get_trajectory.step = step # Save for unit-testing purposes
# Initial spin rate in radians/s
omega0 = 2*pi * rpm0 / 60.
# Glob some constants together
B = -dens*A/(2*m)*h**2
# Set the initial path length to zero
s = 0.
i = 1
while zcur >= get_zground(xcur) and i//step+2 get_zground(xcur) or i//step == 0:
x[i//step+1] = xnext
z[i//step+1] = znext
filled[i//step+1] = True
else: # Don't overwrite the second last point
x[i//step+2] = xnext
z[i//step+2] = znext
filled[i//step+2] = True
# Update the previous and current values for x, z and V
xprev, xcur = xcur, xnext
zprev, zcur = zcur, znext
Vprev, Vcur = Vcur, None
i = i+1 # Increment counter
# Check that the calculation completed
if i//step+2 >= N:
raise RuntimeError, 'Trajectory longer than 15 seconds'
# Return compressed arrays
return x.compress(filled), z.compress(filled)
def get_x_landing(x,z,get_zground=get_zground_default):
"""Determines exactly where the projectile strikes ground.
x,z - the coordinates arrays decribing the trajectory
Returns the x position value.
"""
if x[-1]==x[-2] or get_zground(x[-2])==z[-2]:
# Special cases: 1) coming straight down; and 2) second last point
# is at ground
return x[-2]
else:
# Find z=mx+b between last two points in trajectory
m = (z[-1]-z[-2])/(x[-1]-x[-2])
b = (z[-1]*x[-2]-z[-2]*x[-1]) / (x[-2]-x[-1])
# Determine where the trajectory hits ground
return scipy.optimize.brentq(lambda x: (m*x+b)-get_zground(x),
x[-2], x[-1], xtol=0.1)
def get_roll(x,z,get_zground=get_zground_default,alpha=None):
"""Determines roll coordinates of the golf ball for up to 5 seconds.
x,z - the coordinates arrays decribing the trajectory
get_zground(x) - a function that returns the ground altitude (m) for
position x
alpha - the roll decay constant (/s); if None then a random value will
be used with a mean value of 1
Returns the x- and z-coordinate arrays decribing the roll.
"""
# Get the location of the ground strike
x0 = get_x_landing(x,z,get_zground)
# Initial speed
V0 = (x[-1]-x[-2])/h
# Speed decay constant
if alpha == None:
# Calculate using the gamma distribution with parameters k,theta.
# The mean is k*theta. The standard deviation is sqrt(k)*theta.
# Let the parameters depend on the initial ground speed, with a higher
# standard deviation for higher speed. This will thwart bowling.
Vref = 25. # m/s
k = 25.*(Vref/V0)**2
theta = 1./V0
# Note to self: Keep the mean value k*theta equal to one to be
# consistent with the documentation.
alpha = numpy.random.gamma(25.,0.04) # Mean 1, sd 0.2
# Calculate up to 5 seconds of roll
ts = numpy.arange(0.,5.,h)
# Get the roll positions using analytic equations
xs = V0/alpha*(1.-numpy.exp(-alpha*ts)) + x0
zs = numpy.array([get_zground(v) for v in xs])
# Stop when the roll is slower than 1 m/s
cond = numpy.fabs(xs[1:]-xs[:-1])/h > 1.
cond[:2] = True # Ensures there are at least two points
xs = numpy.compress(cond,xs[:-1])
zs = numpy.compress(cond,zs[:-1])
return xs, zs
def get_carry_distance(x,z,get_zground=get_zground_default):
"""Returns the signed horizonal distance (m) traveled in-flight.
x,z - the coordinates arrays decribing the trajectory
get_zground(x) - a function that returns the ground altitude (m) for
position x
"""
xend = get_x_landing(x,z,get_zground)
return xend-x[0]
def get_roll_distance(x,z):
"""Returns the signed horizontal distance (m) rolled after landing.
x,z - the coordinates arrays decribing the roll
"""
return x[-1]-x[0]