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

Physics Fun: How Many Elephants of Stuff Comes Out of the Tailpipe of Your Car?

We fill up our automobiles at the pump daily. Have you ever thought about how much of the stuff your car uses as fuel escapes from the tailpipe every year?  For a bit of fun, let's measure it in elephants per year.

Before we calculate this, we need to know a few conversion factors to help create our elephant equation:

• Gasoline weighs about 6 pounds (lbs) per gallon
• The weight of an elephant is about 10000 lbs per gallon

Now we can use the handy-dandy dimensional analysis we learned in our science classes to create a rough model that cancels out conversion factors. In the below equation q = Odometer reading in miles/week, and u = Mileage in miles/gallon.

? elephants = (q miles/1 week)*(1 gallon/u miles)*(6 lbs/1 gallon)*(1 elephant/10000 lbs)*(52 weeks/1 year)

This conversion factors simplify to a simple little equation: y =0.0312(q/u), where y is the number of elephants per year.

Let's try our model using U.S. national automobile averages.
q = 231 miles/week, u =23 MPG


y=0.0312(q/u) = 0.0312(231/23) = 0.31 elephants/year

So it turns out Americans have almost a third of an elephant of stuff popping out of their exhaust pipes every year on average. Want to calculate your own car's elephants quickly? Use my Elephant Tailpipe Calculator!
Elephant Tailpipe Calculator 1.0