- #1
CrosisBH
- 27
- 4
- Homework Statement
- Recreate Figure 2.7 from Giordano's Computational Physics 2nd Edition
- Relevant Equations
- The stepping functions: (wind is strictly horizontal)
##x_{i+1} = x_i + v_{x,i} \Delta t##
##v_{x,i+1} = v_{x,i} - \frac{B|\vec{v} - \vec{v}_{wind}|(v_{x,i} - v_{wind})}{m} \Delta t##
##y_{i+1} = y_i + v_{y,i} \Delta t##
##v_{y,i+1} = v_{y,i} - \frac{B |\vec{v} - \vec{v}_{wind}|(v_{y,i})}{m} \Delta t - g\Delta t##
Handy drag coefficient equation
##\frac{B}{m} = 0.0039 + \frac{0.0058}{1+\exp\left(\frac{|\vec{v}| - 35m/s}{5m/s}\right)}##
[CODE lang="python" title="My attempt in python" highlight="57,58,61,62"]import matplotlib.pyplot as plt
import numpy as np
#constants and conditions
initial_velocity = 49.1744 #m/s, book has 110mph
velocity_angle = 35 * np.pi / 180 #coverted to radians because numpy only likes radians
gravity = 9.8 #m/s^2
dt=0.1
windspeed = 4.4704 #m/s book has 10mph
#calculating the initial vector components
position = []
velocity = []
position.append([0,1.5])
position_headwind = position
position_tailwind = position
velocity.append([initial_velocity * np.cos(velocity_angle), initial_velocity * np.sin(velocity_angle)])
velocity_headwind = velocity
velocity_tailwind = velocity
def norm (vector):
return (vector[0]**2 + vector[1]**2)**0.5
def calculate_drag_value(velocity):
#This is B/m value given by the book
b_over_m = 0.0039 + 0.0058/(1+np.exp(velocity - 35)/5)
return b_over_m
def calculate_next_velocity(velocity, windspeed):
#velocity[0] is x-comp of velocity
#velocity[1] is y-comp ...
dragged_velocity = norm([velocity[0] - windspeed, velocity[1]])
newx = velocity[0] - calculate_drag_value(norm(velocity)) * dragged_velocity * (velocity[0] - windspeed) * dt
newy = velocity[1] - calculate_drag_value(norm(velocity)) * dragged_velocity * velocity[1] * dt - gravity * dt
return [newx, newy]
def calculate_next_position(position, velocity):
newx = position[0] + velocity[0] * dt
newy = position[1] + velocity[1] * dt
return [newx, newy]
i=0
#this is the python way to do a do-while loop
#e.g. stop the loop if y<=0, physically meaning the ball hit the grown
while True:
velocity.append(calculate_next_velocity(velocity, 0))
#velocity_tailwind.append(calculate_next_velocity(velocity_tailwind, windspeed))
#velocity_headwind.append(calculate_next_velocity(velocity_headwind, -windspeed))
position.append(calculate_next_position(position, velocity))
#position_tailwind.append(calculate_next_position(position_tailwind, velocity_tailwind))
#position_headwind.append(calculate_next_position(position_headwind, velocity_headwind))
i=i+1
if(position[1]<=0):
break
figure = plt.figure(figsize=(5, 5), dpi=80)
#rearranges my list of vectors into two separate lists of x,y v=[[1,2],[3,4]] becomes x = [1,3] y = [2,4]
x,y = zip(*position)
x_tail, y_tail = zip(*position_tailwind)
x_head, y_head = zip(*position_headwind)
plt.plot(x,y)
plt.plot(x_tail, y_tail)
plt.plot(x_head, y_head)
#axis labels
plt.xlabel("x(m)")
plt.ylabel("y(m)")
#setting the ticks
plt.xticks(np.arange(0,151,50))
#axis range
plt.ylim(0,35)
#title of figure
plt.title("Trajectory of a Batted Baseball")
#show the figure (for testing purposes)
plt.show()[/CODE]
which pretty much yields the same plot with no wind.
If I go ahead and uncomment the lines in the while loop. (line 54). The program breaks off into an infinite loop. I know because I outputted the y position of the no-wind and y is always hanging around the 28-31 range and never goes down (even after 1000 iteration). If I comment the lines again, the output is like 30-40 iterations and it's done.
If I change the wind speed of the original path, I can get the headwind and tailwind curves, so my physics is correct, but my programming seems to be at fault here. (So I apologize if this isn't a question for this specific forum) Does anyone have any insight why it infinite loops if I have 3 calculations in the loop but doesn't if it's one? Thank you!