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.

#!/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"

Fun with Python: Electric Field Using Trapezoid Method

Suppose a line of charge is given that lies on the y-axis from y=-10 cm to y=+10 cm. At any point, the electric field can be calculated by the following function:

E_x=(sigma/(4*pi*epsilon_0))*intregral_-10_10(x/((x^2+y^2)^(3/2)))dy

where sigma is the linear charge density.

This python program solves this equation using the trapezoid method for numerical integration and outputs the data into a .dat file.

#Andrew Samuels
#efield.py
#find the electric field for a given function
#Uses the trapezoid method

#import modules
from math import *

#open a data file
dfile = open('efdata.dat','w')

#make some functions
def efunc(x,y):
    return  x/(x**2.+y**2.)**(3/2.)

def efTrap(f,a,b,step=1):
    A = -(f(x,a) + f(x,b))
    S = float(b - a)
    for n in range(step+1):
        y = a + n*S/step
        A += 2*f(x,y)
    return A*float(b-a)/(2*step)

#constant
k = 10.5/1.113e-4 #10.5microC/4*pi*e_0*1e6

#integrate from x=1 to 100 in steps of 5
x=1
while x <= 100 :
    ef = k * efTrap(efunc,-10,10,10000)
    print "x = " + str(x) + " , eField = " + str(ef)
    dfile.write(str(x) + '\t' + str(ef) + '\n')
    x += 5

#close data file
dfile.close()
print "Done! Created efdata.dat"

Fun with Python: Linear Algebra Solver

This program is similar to my Gaussian Elimination script only this time the matrix is solved using Back Substitution, then checked for accuracy using the Sci-Py module.

#Andrew Samuels
#linear equation solver. Does Gauss-Jordon elimination, then back substituion.

from scipy import *
from scipy import linalg
from pdb import *

#function to compute dot product
def dot(a, b):
    return reduce(lambda sum, p: sum + p[0]*p[1], zip(a,b), 0)

matrix = zeros((3,4))
#aug = zeros((3,1))

matrix[0,0] = 1
matrix[0,1] = 1
matrix[0,2] = 1
matrix[0,3] = 0
matrix[1,0] = 1
matrix[1,1] = -2
matrix[1,2] = 2
matrix[1,3] = 4
matrix[2,0] = 1
matrix[2,1] = 2
matrix[2,2] = -1
matrix[2,3] = -1


pivot = 0
ncols = 4
nrows = 3


#for all pivots
for col in range(ncols-1):
    #divide pivot by itself to make it 1
    matrix[pivot,:] = matrix[pivot,:]/matrix[pivot,pivot]
   
    #for all pivots that haven't been done yet
    for row in range(pivot+1,nrows):
        matrix[row,:] = matrix[row,:] - matrix[pivot,:]
    pivot = pivot + 1

print " "
print "Row Echelon Form: "   
print matrix 
print "\n"

#solve the matrix and print solution
for pivot in range(nrows-1,-1,-1):
    matrix[pivot] = (matrix[pivot] - dot(matrix[pivot,pivot+1:nrows],matrix[pivot+1:nrows]))/matrix[pivot,pivot]
print "Solved by Back Subst.: "
print matrix 
print "\n"  

#check with scipy
a = array([[1,1,1],[1,-2,2],[1,2,-1]])
b = array([[0],[4],[-1]])
sol = linalg.solve(a,b)
print "Solved by Scipy (check): "
print sol
print "\n"

Fun with Python: Gaussian Elimination

Here is my python program that does Gaussian Elimination on a 3 x 4 matrix. Of course, there are existing modules with functions that could do this automatically, but this was a good learning exercise doing it the hard way.

#Andrew Samuels
#Gaussian Elimination
#performs gaussian elimination on a 3 by 4 matrix

#import scipy module
from scipy import *
from pdb import *

#construct the matrix
matrix = zeros((3,4))

matrix[0,0] = 1
matrix[0,1] = 1
matrix[0,2] = 1
matrix[0,3] = 0
matrix[1,0] = 1
matrix[1,1] = -2
matrix[1,2] = 2
matrix[1,3] = 4
matrix[2,0] = 1
matrix[2,1] = 2
matrix[2,2] = -1
matrix[2,3] = -1


pivot = 0
ncols = 4
nrows = 3

#for all pivots
for col in range(ncols-1):
    #divide pivot by itself to make it 1
    matrix[pivot,:] = matrix[pivot,:]/matrix[pivot,pivot]
   
    #for all pivots that haven't been done yet
    for row in range(pivot+1,nrows):
        matrix[row,:] = matrix[row,:] - matrix[pivot,:]
    pivot = pivot + 1
   

print matrix

Fun with Python: Factorial Calculator

This program written in python calculates a factorial when given an integer. The program will complain if you don't give it a valid integer, such as a letter or dollar sign.

#Andrew Samuels
#factorial calculator
#Fall 2011
#calculates a factorial from integer input from user

n = raw_input('Enter an integer: ')
try:
    n=int(n)
except:
    print 'Enter a valid integer, dummy!'
    exit(1)

f=1
while n > 1 :
    f = f * n
    n = n - 1
    print f