#!/opt/local/bin/python
#Andrew Samuels
#baseball.py
from scipy import *
from pdb import *
#create arrays
dt = 1e-4
nMAX = int(1/dt)
x = zeros(nMAX)
y = zeros(nMAX)
z = zeros(nMAX)
vx = zeros(nMAX)
vy = zeros(nMAX)
vz = zeros(nMAX)
v = zeros(nMAX)
#constants
B = 4.1e-4 #magnus force (dimensionless)
w = 1800/60.0 * 2*pi #rotaton rate of baseball, 1800 rpm to rad/s
g = 9.81 #acceleration of gravity (m/s^2)
v0_fb = 42.4688 #initial speed of pitch, fastball (m/s)
v0_op = 38.0 #initial speed of pitch, other pitches (m/s)
l = 18.44 #distance from pitcher (m)
T_fb = l/v0_fb #step-length, estimated flight time, fastball
T_op = l/v0_op #step-length, estimated flight time, other pitches
nT_fb = int(T_fb/dt) #number of iterations for fastball
nT_op = int(T_op/dt) #number of iterations for other pitches
d2r = pi/180 # convert degrees to radians
theta = 1 * d2r # elevation angle of pitch, degrees
phi1 = 225 * d2r #direction of w, relative to z, fastball
phi2 = 45 * d2r #direction of w, relative to z, curveball
phi3 = 0 * d2r #direction of w, relative to z, slider
phi4 = 135 * d2r #direction of w, relative to z, screwball
#drag force constants
c1 = 0.0039
c2 = 0.0058
vd = 35.0 #(m/s)
delta = 5.0 # (m/s)
#initial conditions
x[0] = 0.0
y[0] = 0.0
z[0] = 0.0
vx[0] = v0_fb * cos(theta)
vy[0] = 0.0
vz[0] = v0_fb * sin(theta)
#pass dx/dt=vx through runkut4
def rk4dx(vx,dt):
def dxdt(vx):
dxdt = vx
return dxdt
K0 = dt * dxdt(vx)
K1 = dt * dxdt(vx + K0/2.0)
K2 = dt * dxdt(vx + K1/2.0)
K3 = dt * dxdt(vx + K2)
vxnew = (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
return vxnew
#pass dy/dt=vy through runkut4
def rk4dy(vy,dt):
def dydt(vy):
dydt = vy
return dydt
K0 = dt * dydt(vy)
K1 = dt * dydt(vy + K0/2.0)
K2 = dt * dydt(vy + K1/2.0)
K3 = dt * dydt(vy + K2)
vynew = (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
return vynew
#pass dz/dt=vz through runkut4
def rk4dz(vz,dt):
def dzdt(vz):
dzdt = vz
return dzdt
K0 = dt * dzdt(vz)
K1 = dt * dzdt(vz + K0/2.0)
K2 = dt * dzdt(vz + K1/2.0)
K3 = dt * dzdt(vz + K2)
vznew = (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
return vznew
#solve accelerations with runge kutta 4
def rk4(v,vx,vy,vz,dt,phi):
def ax(v,vx,vy,vz,phi):
#drag force function, Fd(v)
def Fd(v):
return c1 + (c2/(1.0 + exp((v-vd)/delta)))
return -Fd(v)*v*vx + B*w*(vz*sin(phi)-vy*cos(phi))
def ay(v,vx,vy,phi):
#drag force function, Fd(v)
def Fd(v):
return c1 + (c2/(1.0 + exp((v-vd)/delta)))
return -Fd(v)*v*vy + B*w*vx*cos(phi1)
def az(v,vx,vz,phi):
#drag force function, Fd(v)
def Fd(v):
return c1 + (c2/(1.0 + exp((v-vd)/delta)))
return -g*-Fd(v)*v*vz - B*w*vx*sin(phi)
K0 = dt * ax(v,vx,vy,vz,phi)
L0 = dt * ay(v,vx,vy,phi)
M0 = dt * az(v,vx,vz,phi)
K1 = dt * ax(v + dt/2.0,vx + K0/2.0,vy + L0/2.0,vz + M0/2.0,phi)
L1 = dt * ay(v + dt/2.0,vx + K0/2.0,vy + L0/2.0,phi)
M1 = dt * az(v + dt/2.0,vx + K0/2.0,vz + M0/2.0,phi)
K2 = dt * ax(v + dt/2.0,vx + K1/2.0,vy + L1/2.0,vz + M1/2.0,phi)
L2 = dt * ay(v + dt/2.0,vx + K1/2.0,vy + L1/2.0,phi)
M2 = dt * az(v + dt/2.0,vx + K1/2.0,vz + M1/2.0,phi)
K3 = dt * ax(v + dt,vx + K2,vy + L2,vz + M0,phi)
L3 = dt * ay(v + dt,vx + K1,vy + L2,phi)
M3 = dt * az(v + dt,vx + K1,vz + M2,phi)
axnew = (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
aynew = (L0 + 2.0*L1 + 2.0*L2 + L3)/6.0
aznew = (M0 + 2.0*M1 + 2.0*M2 + M3)/6.0
return axnew,aynew,aznew
f=open('fastball.dat','w')
#first time step is n + 1
#start loop
for i in range(0,nT_fb):
v = (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i])
v = (v)**0.5
print "X: " + str(x[i]) + " , Y: " + str(y[i]) + " , Z: " + str(z[i]) + " ,T: " + str(x[i]/v) +" , V: "+ str(v) + "\n"
f.write(str(x[i])+'\t'+str(y[i])+'\t'+ str(z[i])+'\t'+str(x[i]/v) + '\t' + str(v) +'\n')
x[i+1] = x[i] + rk4dx(vx[i],dt)
y[i+1] = y[i] + rk4dy(vy[i],dt)
z[i+1] = z[i] + rk4dz(vz[i],dt)
vx[i+1] = vx[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi1)[0]
vy[i+1] = vy[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi1)[1]
vz[i+1] = vz[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi1)[2]
f.close()
print "Created fastball.dat"
#change initial conditions for the other pitches
vx[0] = v0_op * cos(theta)
vz[0] = v0_op * sin(theta)
#curveball
h=open('curveball.dat','w')
for i in range(0,nT_op):
v = (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i])
v = (v)**0.5
print "X: " + str(x[i]) + " , Y: " + str(y[i]) + " , Z: " + str(z[i]) + " ,T: " + str(x[i]/v) +" , V: "+ str(v) + "\n"
h.write(str(x[i])+'\t'+str(y[i])+'\t'+ str(z[i])+'\t'+str(x[i]/v) + '\t' + str(v) +'\n')
x[i+1] = x[i] + rk4dx(vx[i],dt)
y[i+1] = y[i] + rk4dy(vy[i],dt)
z[i+1] = z[i] + rk4dz(vz[i],dt)
vx[i+1] = vx[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi2)[0]
vy[i+1] = vy[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi2)[1]
vz[i+1] = vz[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi2)[2]
h.close()
print "Created curveball.dat"
#slider
j=open('slider.dat','w')
for i in range(0,nT_op):
v = (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i])
v = (v)**0.5
print "X: " + str(x[i]) + " , Y: " + str(y[i]) + " , Z: " + str(z[i]) + " ,T: " + str(x[i]/v) +" , V: "+ str(v) + "\n"
j.write(str(x[i])+'\t'+str(y[i])+'\t'+ str(z[i])+'\t'+str(x[i]/v) + '\t' + str(v) +'\n')
x[i+1] = x[i] + rk4dx(vx[i],dt)
y[i+1] = y[i] + rk4dy(vy[i],dt)
z[i+1] = z[i] + rk4dz(vz[i],dt)
vx[i+1] = vx[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi3)[0]
vy[i+1] = vy[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi3)[1]
vz[i+1] = vz[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi3)[2]
j.close()
print "Created slider.dat"
#screwball
o=open('screwball.dat','w')
for i in range(0,nT_op):
v = (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i])
v = (v)**0.5
print "X: " + str(x[i]) + " , Y: " + str(y[i]) + " , Z: " + str(z[i]) + " ,T: " + str(x[i]/v) +" , V: "+ str(v) + "\n"
o.write(str(x[i])+'\t'+str(y[i])+'\t'+ str(z[i])+'\t'+str(x[i]/v) + '\t' + str(v) +'\n')
x[i+1] = x[i] + rk4dx(vx[i],dt)
y[i+1] = y[i] + rk4dy(vy[i],dt)
z[i+1] = z[i] + rk4dz(vz[i],dt)
vx[i+1] = vx[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi4)[0]
vy[i+1] = vy[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi4)[1]
vz[i+1] = vz[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi4)[2]
o.close()
print "Created screwball.dat"
20111228
Fun with Python: Simulating Baseball Pitches using 4th Order Runge-Kutta Method
Here is a description of the problem in PDF format. Basically this code simulates four different baseball pitches (slider, curve, fastball, screwball) by solving differential equations using the 4th order Runge-Kutta method.
