在 python 中修复射弹运动公式

Fixing projectile motion formulas in python

该代码输出的图形与您对抛射运动类型图形的期望相去甚远。此外,如果将步数更改为 2 左右,则图表不会输出太多内容。

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

看起来你没有更新你的角度,所以它没有达到零 vx,而是向后加速,因为它认为你仍然以 60 度的初始角度前进。添加基于当前速度矢量的更新角度结果如下:

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)")