# Help in Python for Computational Physics

Hi, I am in a computational physics class and we are programming projectile motion in python. I do not know why I am having such a hard time at this, but I cannot get my programs to work.

1. Homework Statement
Calculate the effect of backspin on a fastball. How much does an angular velocity of 1000 rpm affect the trajectory. We are incorporating both drag and spin into the equations. I am the Euler method to solve the problem.

2. Homework Equations
Drag: F(v) = 0.0039 + 0.0058/(1 + exp((v-35)/5)) (As an aside, there is a graph in my book modeling this equation and I cannot replicate it. It shows that at about 25 mph the function starts at .5 and at 55 or so, there is an exponential drop to .2)

dx/dt = vx
dy/dt = vy ,
dz/dt = vz ,

dvx/dt = −F(v) v vx + B ω (vz sin φ − vy cos φ),

dvy/dt = −F(v) v vy + B ω vx cos φ,

dvz/dt = −g − F(v) v vz − B ω vx sin φ.

v = 44.7 m/s
B = 4.1*10^-4 (unitless)
φ = 225 degrees

x(t = 0) = 0,
y(t = 0) = 0,
z(t = 0) = 3,
vx (t = 0) = v0 cos θ
vy (t = 0) = 0,
vz (t = 0) = v0 sin θ,

3. The Attempt at a Solution

Below is the code that I have so far.

from scipy import *
from pylab import *
from numpy import *

#Define Constants
m = .149 #mass of baseball in kilograms
B = (4.1*10**-4) #Spin coefficient - unitless
w = (100*pi)/3 #1000 rpm in rad/s
g = 9.81 #acceleration due to gravity m/s^2

theta = 1
pho = 225

#Define Arrays
dt = 1*10**-4
n = int(1/dt)
phi = pi/180 #Convert degrees to radians

t = linspace(0,0,n+1)
x = linspace(0,0,n+1)
y = linspace(0,0,n+1)
z = linspace(0,0,n+1)
vx = linspace(0,0,n+1)
vy = linspace(0,0,n+1)
vz = linspace(0,0,n+1)
v = linspace(0,0,n+1)
F = linspace(0,0,n+1)

#Initial Conditions
x = 0 #Axis is oriented os that x is the line from the pitcher's mound to homeplate
y = 0
z = 3
v = 44.704 #Initial speed of the ball in m/s (100mph)
F = 0.0039 + (0.0058)/(1 + exp((v - 35)/5))

vx = v*cos(theta*phi) #the ball is initially thrown at a small 1 degree angle above the horizontal
vy = 0 #No inital velocity in the y direction
vz = v*sin(theta*phi)
i = 0
for i in xrange(n):
t[i+1] = t + dt
v[i+1] = sqrt(vx**2 + vy**2 + vz**2)
F[i+1] = 0.0039 + (0.0058)/(1 + exp((v - 35)/5))
vx[i+1] = -F * v * vx + B * w *(vz*sin(pho*phi) - vy*cos(pho*phi))
x[i+1] = x + vx*dt
y[i+1] = y + vy*dt
vy[i+1] = -F * v * vy + B * w * vx * cos(pho*phi)
z[i+1] = z + vz*dt
vz[i+1] = -g - F * v * vz - B * w * vx * sin(pho*phi)
print(v, x, y, z, F, t)
if z < 0:
break