Fixing projectile motion formulas in python

Question

The code outputs a graph that is nowhere near what you'd expect of a projectile motion type of graph. In addition, if you change the step number to around 2, the graph doesn't output much of anything.

import numpy as np
import matplotlib.pyplot as plt

grav = 9.8
airDen = 1.225      # kg/m^3
areaProj = 0.8643   # m^2
v0 = 198.1          # m/s
angle = 60.0   
dragCoef = 0.62     #Approximate value
mass = 15.0         #kg
step = 0.02    #time step in seconds 

t = [0]
vx = [v0*np.cos(angle*np.pi/180)]
vy = [v0*np.sin(angle*np.pi/180)]
x = [0]
y = [0]

dragForce = 0.5*dragCoef*areaProj*(v0**2)*airDen


ax = [-(dragForce*np.cos(angle/180*np.pi))/mass]
ay = [-grav-(dragForce*np.sin(angle/180*np.pi)/mass)]


counter = 0
while(y[counter] >= 0):
    t.append(t[counter]+step)

    vx.append(vx[counter]+step*ax[counter])
    vy.append(vy[counter]+step*ay[counter])

    x.append(x[counter]+step*vx[counter])
    y.append(y[counter]+step*vy[counter])

    vel = np.sqrt(vx[counter+1]**2 + vy[counter+1]**2)
    dragForce = 0.5*dragCoef*areaProj*(vel**2)*airDen

    ax.append(-(dragForce*np.cos(angle/180*np.pi))/mass)
    ay.append(-grav-(dragForce*np.sin(angle/180*np.pi)/mass))

    counter=counter+1

plt.plot(x,y)
plt.ylabel("y (m)")
plt.xlabel("x (m)")

Graph-example

score 2 · Accepted Answer · answered Dec 09 '19 at 23:56

Looks like you're not updating your angle, so instead of reaching zero vx it accelerates backwards since it thinks you're still going at 60 degrees initial angle. Adding an updated angle based on the current velocity vector results in the following:

import numpy as np
import matplotlib.pyplot as plt

grav = 9.8
airDen = 1.225      # kg/m^3
areaProj = 0.8643   # m^2
v0 = 198.1          # m/s
angle = 60.0   
dragCoef = 0.62     #Approximate value
mass = 15.0         #kg
step = 0.02    #time step in seconds 

t = [0]
vx = [v0*np.cos(angle*np.pi/180)]
vy = [v0*np.sin(angle*np.pi/180)]
x = [0]
y = [0]

dragForce = 0.5*dragCoef*areaProj*(v0**2)*airDen


ax = [-(dragForce*np.cos(angle/180*np.pi))/mass]
ay = [-grav-(dragForce*np.sin(angle/180*np.pi)/mass)]


counter = 0
while(y[counter] >= 0):
    t.append(t[counter]+step)

    vx.append(vx[counter]+step*ax[counter])
    vy.append(vy[counter]+step*ay[counter])

    x.append(x[counter]+step*vx[counter])
    y.append(y[counter]+step*vy[counter])

    vel = np.sqrt(vx[counter+1]**2 + vy[counter+1]**2)
    dragForce = 0.5*dragCoef*areaProj*(vel**2)*airDen

    angle = np.arctan(vy[counter]/vx[counter]) * (180 / np.pi)

    ax.append(-(dragForce*np.cos(angle/180*np.pi))/mass)
    ay.append(-grav-(dragForce*np.sin(angle/180*np.pi)/mass))

    counter=counter+1


plt.plot(x,y)
plt.ylabel("y (m)")
plt.xlabel("x (m)")

Thank you! Would you happen to know why changing the step to 1 or close to 1, alters the graph? — hero7209, Dec 10 '19 at 01:02
Also, given those numbers I should be seeing a range of around 3,500 m on x. This number is only achievable when I take out the velocity factor out of the drag force calculation. — hero7209, Dec 10 '19 at 01:20
Changing the timestep makes the simulation more coarse, and for something like this you should be using an ODE integrator or your own integration scheme. Check out the offerings from SciPy here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html — Philip Ciunkiewicz, Dec 10 '19 at 01:57
Also, if the answer resolved your question check it off as the accepted answer :) — Philip Ciunkiewicz, Dec 10 '19 at 02:03

Fixing projectile motion formulas in python

1 Answers1