baseball_trajectory.py

import math
import matplotlib.pyplot as plt

#initial condition
mass = 5.125 #Oz
circumference = 9.125 #inches
x0  = 0 #ft
y0 = 2.000 #ft
z0 = 3 #ft
exit_speed  = 100 #mph
launch_angle = 30 #degree
direction  = 0 #degree
sign = 1
backspin= -763+120*launch_angle+21*direction*sign#rpm
sidespin = -sign*849-94*direction #rpm
wg = 0 #rpm
tau = 25 #sec
dt = 0.01


Temp_F = 78 #F
elev_ft = 0 #ft
vwind = 0 #mph
phiwind  = 0 #degree
hwind = 0 #mph
relative_humidity = 50 #%
barometric_pressure = 29.92 #Hg

Temp_C = (5/9)*(Temp_F-32)
elev_m = elev_ft / 3.2808

vxw = vwind*1.467*math.sin(phiwind*math.pi/180) #ft/s
vyw = vwind*1.467*math.cos(phiwind*math.pi/180) #ft/s

beta = 0.0001217 #1/m
SVP = 4.5841*math.exp((18.687-Temp_C /234.5)*Temp_C/(257.14+Temp_C ))
barometric_pressure_mmHg = barometric_pressure*1000/39.37
rho_kg_m3 =  1.2929*(273/(Temp_C+273)*(barometric_pressure_mmHg*math.exp(-beta*elev_ft)-0.3783*relative_humidity*SVP/100)/760) #kg/m^3
rho_lb_ft3 = rho_kg_m3*0.06261 #lb/ft^3
# const
c0 = 0.07182*rho_lb_ft3*(5.125/mass)*(circumference/9.125)**2

v0 = 1.467*exit_speed
v0x = v0*math.cos(launch_angle*math.pi/180)*math.sin(direction*math.pi/180)
v0y = v0*math.cos(launch_angle*math.pi/180)*math.cos(direction*math.pi/180)
v0z = v0*math.sin(launch_angle*math.pi/180)
wx = (backspin*math.cos(direction*math.pi/180)- \
      sidespin*math.sin(launch_angle*math.pi/180)*math.sin(direction*math.pi/180)+wg*v0x/v0)*math.pi/30
wy = (-backspin*math.sin(direction*math.pi/180)- \
      sidespin*math.sin(launch_angle*math.pi/180)*math.cos(direction*math.pi/180)+wg*v0y/v0)*math.pi/30
wz = (sidespin*math.cos(launch_angle*math.pi/180)+wg*v0z/v0)*math.pi/30
omega = math.sqrt(wx*wx+wy*wy+wz*wz)
romega = (circumference/2/math.pi)*omega/12



def accelaration(t,vx,vy,vz,z):
    v = math.sqrt(vx*vx+vy*vy+vz*vz)
    if(z > hwind):
        vw = math.sqrt((vx-vxw)**2 + (vy-vyw)**2+vz**2)
        vyw_sim = vyw
    else:
        vw = v
        vyw_sim = 0
    S = (romega/vw)*math.exp(-t/(tau*146.7/v))    
    Cd = 0.4105*(1+0.2017*S*S)
    Cl = 1/(2.32+0.4/S)
    vxw_sim = 0
    adragx = -c0*Cd*vw*(vx-vxw)
    adragy = -c0*Cd*vw*(vy-vyw_sim)
    adragz = -c0*Cd*vw*vz
    w = omega*math.exp(-t/tau)*30/math.pi
    aMagx = c0*(Cl/omega)*vw*(wy*vz-wz*(vy-vyw))
    aMagy = c0*(Cl/omega)*vw*(wz*(vx-vw)-wx*vz)
    aMagz = c0*(Cl/omega)*vw*(wx*(vy-vyw)-wy*(vx-vxw))
    ax = adragx+aMagx
    ay = adragy+aMagy
    az = adragz+aMagz-32.174

    
    return ax,ay,az

def RK4(x,y,z,vx,vy,vz,t):
    k1x,k1y,k1z = tuple([dt*i for i in accelaration(t,vx,vy,vz,z)]) #dt * f(x,t)
    k2x,k2y,k2z = tuple([dt*i for i in accelaration(t+dt/2,vx+k1x/2,vy+k1y/2,vz+k1z/2,z)])
    k3x,k3y,k3z = tuple([dt*i for i in accelaration(t+dt/2,vx+k2x/2,vy+k2y/2,vz+k2z/2,z)])
    k4x,k4y,k4z = tuple([dt*i for i in accelaration(t+dt,vx+k3x,vy+k3y,vz+k3z,z)])
       
    
    vx += (1/6)*(k1x+2*k2x+2*k3x+k4x)
    vy += (1/6)*(k1y+2*k2y+2*k3y+k4y)
    vz += (1/6)*(k1z+2*k2z+2*k3z+k4z)
    
    x += vx * dt
    y += vy * dt
    z += vz * dt

    return x,y,z,vx,vy,vz

# y_path_euler = []
# z_path_euler = []

# x = 0
# y = y0
# z = z0

# vx = v0x
# vy = v0y
# vz = v0z
# for t in range(1000):
    
#     y_path_euler.append(y)
#     z_path_euler.append(z)
    
#     t = dt*t

#     ax, ay, az = accelaration(t,vx,vy,vz,z)
    
#     y += vy * dt + ay * dt * dt
#     z +=  vz * dt + az * dt * dt

#     vy += ay * dt
#     vz += az * dt
    
    
#     if(z <0):

#         print('Hang time = %1.2f s Distance = %1.1f feet' % (t, y))
        
#         break
x_path_RK4 = []
y_path_RK4 = []
z_path_RK4 = []

x = 0
y = y0
z = z0

vx = v0x
vy = v0y
vz = v0z
for t in range(1000):
    
    x_path_RK4.append(x)
    y_path_RK4.append(y)
    z_path_RK4.append(z)
    
    t = dt*t
    
    x,y,z,vx,vy,vz = RK4(x,y,z,vx,vy,vz,t)
    
    if(z <0):

        print('Hang time = %1.2f s Distance = %1.1f feet' % (t, y))
        
        break     
            
            
#%matplotlib inline

# plt.plot(y_path_RK4,z_path_RK4, label='RK4')
# #plt.plot(y_path_euler,z_path_euler, label='Euler')
# #plt.gca().legend(('Euler','RK4'))
# plt.title('baseball tracjectory')
# plt.xlabel('horizontal distance (ft)')
# plt.ylabel('vertical distance (ft)')
# plt.grid(True)
# plt.axes().set_aspect('equal', 'datalim') 

from mpl_toolkits.mplot3d import axes3d, Axes3D

fig = plt.figure()
ax = Axes3D(fig)
ax.plot3D(x_path_RK4,y_path_RK4,z_path_RK4)
plt.xlim([-200,200])
plt.ylim([0,400])
plt.title('baseball tracjectory')