我的代码挂在哪里（无限循环中）？

发布于 2025-01-10 20:12:31 字数 6123 浏览 0 评论 0原文

我是 Python 新手，试图让这个脚本运行，但它似乎挂在无限循环中。当我使用 ctrl+c 停止它时，它总是在第 103 行。

vs = 20.05 * np.sqrt(Tb + Lb * (y - y0)) # m/s speed of sound as a function of temperature

我习惯了 MatLab（来自学校）及其编辑器。我之前遇到过此代码的编码问题。对（免费）编辑器有什么建议吗？我目前正在使用 JEdit 和/或记事本。

这是完整的脚本：

#!/usr/bin/env python
# -*- coding: ANSI -*-
import numpy as np 
from math import *
from astropy.table import Table 
import matplotlib.pyplot as plt
from hanging_threads import start_monitoring#test for code hanging
start_monitoring(seconds_frozen=10, test_interval=100)

"""Initial Conditions and Inputs"""

d = 154.71/1000 # diameter of bullet (in meters) 
m = 46.7 # mass of bullet ( in kg)
K3 = 0.87*0.3735 # drag coefficient at supersonic speed 
Cd1 = 0.87*0.108 #drag coefficient at subsonic speed
v0 = 802 # muzzle velocity in m/sec 
dt = 0.01 # timestep in seconds

"""coriolis inputs"""
L = 90*np.pi/180 # radians - latitude of firing site
AZ = 90*np.pi/180 # radians - azimuth angle of fire measured clockwise from North
omega = 0.0000727 #rad/s rotation of the earth

"""wind inputs""" 
wx = 0 # m/s
wz = 0  # m/s

"""initializing variables""" 
vx = 0 #initial x velocity 
vy = 0 #initial y velocity 
vy0 = 0
y_max = 0 #apogee 
v = 0
t = 0
x = 0

"""Variable Atmospheric Pressure"""
rho0 = 1.2041 # density of air at sea-level (kg/m^3) 
T = 20 #temperature at sea level in celcius
Tb = T + 273.15 # temperature at sea level in Kelvin
Lb = -2/304.8 # temperature lapse rate in K/m (-2degrees/1000ft)- not valid above 36000ft
y = 0 # current altitude 
y0 = 0 # initial altitude
g = 9.81 # acceleration due to gravity in m/s/s 
M = 0.0289644 #kg/mol # molar mass of air
R = 8.3144598 # J/molK - universal gas constant

# air density as a function of altitude and temperature 
rho = rho0 * ((Tb/(Tb+Lb*(y-y0)))**(1+(g*M/(R*Lb))))

"""Variable Speed of Sound"""

vs = 20.05*np.sqrt(Tb +Lb*(y-y0)) # m/s speed of sound as a function of temperature

Area = pi*(d/2)**2 # computing the reference area 
phi_incr = 5 #phi0 increment (degrees)
N = 12  # length of table

"""Range table"""
dtype = [('phi0', 'f8'), ('phi_impact', 'f8'), ('x', 'f8'), ('z', 'f8'),('y', 'f8'), ('vx', 'f8'), ('vz', 'f8'), ('vy', 'f8'), ('v', 'f8'),('M', 'f8'), ('t', 'f8')]
table = Table(data=np.zeros(N, dtype=dtype))

"""Calculates entire trajectory for each specified angle""" 
for i in range(N):
    phi0 = (i + 1) * phi_incr

    """list of initial variables used in while loop""" 
    t = 0
    y = 0
    y_max = y
    x = 0
    z = 0
    vx = v0*np.cos(radians(phi0))
    vy = v0*np.sin(radians(phi0))
    vx_w = 0
    vz_w = 0
    vz = 0
    v = v0
    ay = 0
    ax = 0
    wx = wx
    wz = wz
    rho = rho0 * ((Tb / (Tb + Lb * (y - y0))) ** (1 + (g * M / (R * Lb))))
    vs = 20.05 * np.sqrt(Tb + Lb * (y - y0)) # m/s speed of sound as a function of temperature
    ax_c = -2 * omega * ((vz * sin(L)) + vy * cos(L) * sin(AZ))
    ay_c = 2 * omega * ((vz * cos(L) * cos(AZ)) + vx_w * cos(L) * sin(AZ))
    az_c = -2 * omega * ((vy * cos(L) * cos(AZ)) - vx_w * sin(L))
    Mach = v/vs
    
    """ initializing variables for plots"""
    t_list = [t]
    x_list = [x]
    y_list = [y]
    vy_list = [vy]
    v_list = [v]
    phi0_list = [phi0]
    Mach_list = [Mach]
    
    while y >= 0:
        phi0 = phi0
        """drag calculation with variable density, Temp and sound speed"""
        rho = rho0 * ((Tb / (Tb + Lb * (y - y0))) ** (1 + (g * M / (R *Lb))))
        vs = 20.05 * np.sqrt(Tb + Lb * (y - y0)) # m/s speed of sound as a function of temperature
        Cd3 = K3 / sqrt(v / vs)
        Mach = v/vs
        
        """Determining drag regime"""
        if v > 1.2 * vs: #supersonic
            Cd = Cd3
        elif v < 0.8 * vs: #subsonic
            Cd = Cd1
        else: #transonic
            Cd = ((Cd3 - Cd1)*(v/vs - 0.8)/(0.4)) + Cd1
            
        """Acceleration due to Coriolis"""
        ax_c = -2*omega*((vz_w*sin(L))+ vy*cos(L)*sin(AZ))
        ay_c = 2*omega*((vz_w*cos(L)*cos(AZ))+ vx_w*cos(L)*sin(AZ))
        az_c = -2*omega*((vy*cos(L)*cos(AZ))- vx_w*sin(L))

        """Total acceleration calcs"""
        if vx > 0:
            ax = -0.5*rho*((vx-wx)**2)*Cd*Area/m + ax_c
        else:
            ax = 0

        """ Vy before and after peak"""
        if vy > 0:
            ay = (-0.5 * rho * (vy ** 2) * Cd * Area / m) - g + ay_c
        else:
            ay = (0.5 * rho * (vy ** 2) * Cd * Area / m) - g + ay_c
            az = az_c
            vx = vx + ax*dt # vx without wind
            # vx_w = vx with drag and no wind + wind
            vx_w = vx + 2*wx*(1-(vx/v0*np.cos(radians(phi0))))
            vy = vy + ay*dt
            vz = vz + az*dt
            vz_w = vz + wz*(1-(vx/v0*np.cos(radians(phi0))))
            """projectile velocity"""
            v = sqrt(vx_w**2 + vy**2 + vz**2)
            """new x, y, z positions"""
            x = x + vx_w*dt
            y = y + vy*dt
            z = z + vz_w*dt
            if y_max <= y:
                y_max = y
                phi_impact = degrees(atan(vy/vx)) #impact angle in degrees
                """ appends selected data for ability to plot"""
                t_list.append(t)
                x_list.append(x)
                y_list.append(y)
                vy_list.append(vy)
                v_list.append(v)
                phi0_list.append(phi0)
                Mach_list.append(Mach)
                if y < 0:
                    break
                    t += dt
                    
            """Range table output"""
            table[i] = ('%.f' % phi0, '%.3f' % phi_impact, '%.1f' % x,'%.2f' % z, '%.1f' % y_max, '%.1f' % vx_w,'%.1f' % vz,'%.1f' % vy,'%.1f' % v,'%.2f' %Mach, '%.1f' % t)

    """ Plot"""
    plt.plot(x_list, y_list, label='%d°' % phi0)#plt.plot(x_list, y_list, label='%d°' % phi0)
    plt.title('Altitude versus Range')
    plt.ylabel('Altitude (m)')
    plt.xlabel('Range (m)')
    plt.axis([0, 30000, 0, 15000])
    plt.grid(True)
    
print(table)

legend = plt.legend(title="Firing Angle",loc=0, fontsize='small', fancybox=True)

plt.show()

提前谢谢您

原文

I am new to Python and trying to get this script to run, but it seems to be hanging in an infinite loop. When I use ctrl+c to stop it, it is always on line 103.

vs = 20.05 * np.sqrt(Tb + Lb * (y - y0)) # m/s speed of sound as a function of temperature

I am used to MatLab (from school) and the editor it has. I ran into issues earlier with the encoding for this code. Any suggestions on a (free) editor? I am currently using JEdit and/or Notepad.

Here is the full script:

#!/usr/bin/env python
# -*- coding: ANSI -*-
import numpy as np 
from math import *
from astropy.table import Table 
import matplotlib.pyplot as plt
from hanging_threads import start_monitoring#test for code hanging
start_monitoring(seconds_frozen=10, test_interval=100)

"""Initial Conditions and Inputs"""

d = 154.71/1000 # diameter of bullet (in meters) 
m = 46.7 # mass of bullet ( in kg)
K3 = 0.87*0.3735 # drag coefficient at supersonic speed 
Cd1 = 0.87*0.108 #drag coefficient at subsonic speed
v0 = 802 # muzzle velocity in m/sec 
dt = 0.01 # timestep in seconds

"""coriolis inputs"""
L = 90*np.pi/180 # radians - latitude of firing site
AZ = 90*np.pi/180 # radians - azimuth angle of fire measured clockwise from North
omega = 0.0000727 #rad/s rotation of the earth

"""wind inputs""" 
wx = 0 # m/s
wz = 0  # m/s

"""initializing variables""" 
vx = 0 #initial x velocity 
vy = 0 #initial y velocity 
vy0 = 0
y_max = 0 #apogee 
v = 0
t = 0
x = 0

"""Variable Atmospheric Pressure"""
rho0 = 1.2041 # density of air at sea-level (kg/m^3) 
T = 20 #temperature at sea level in celcius
Tb = T + 273.15 # temperature at sea level in Kelvin
Lb = -2/304.8 # temperature lapse rate in K/m (-2degrees/1000ft)- not valid above 36000ft
y = 0 # current altitude 
y0 = 0 # initial altitude
g = 9.81 # acceleration due to gravity in m/s/s 
M = 0.0289644 #kg/mol # molar mass of air
R = 8.3144598 # J/molK - universal gas constant

# air density as a function of altitude and temperature 
rho = rho0 * ((Tb/(Tb+Lb*(y-y0)))**(1+(g*M/(R*Lb))))

"""Variable Speed of Sound"""

vs = 20.05*np.sqrt(Tb +Lb*(y-y0)) # m/s speed of sound as a function of temperature

Area = pi*(d/2)**2 # computing the reference area 
phi_incr = 5 #phi0 increment (degrees)
N = 12  # length of table

"""Range table"""
dtype = [('phi0', 'f8'), ('phi_impact', 'f8'), ('x', 'f8'), ('z', 'f8'),('y', 'f8'), ('vx', 'f8'), ('vz', 'f8'), ('vy', 'f8'), ('v', 'f8'),('M', 'f8'), ('t', 'f8')]
table = Table(data=np.zeros(N, dtype=dtype))

"""Calculates entire trajectory for each specified angle""" 
for i in range(N):
    phi0 = (i + 1) * phi_incr

    """list of initial variables used in while loop""" 
    t = 0
    y = 0
    y_max = y
    x = 0
    z = 0
    vx = v0*np.cos(radians(phi0))
    vy = v0*np.sin(radians(phi0))
    vx_w = 0
    vz_w = 0
    vz = 0
    v = v0
    ay = 0
    ax = 0
    wx = wx
    wz = wz
    rho = rho0 * ((Tb / (Tb + Lb * (y - y0))) ** (1 + (g * M / (R * Lb))))
    vs = 20.05 * np.sqrt(Tb + Lb * (y - y0)) # m/s speed of sound as a function of temperature
    ax_c = -2 * omega * ((vz * sin(L)) + vy * cos(L) * sin(AZ))
    ay_c = 2 * omega * ((vz * cos(L) * cos(AZ)) + vx_w * cos(L) * sin(AZ))
    az_c = -2 * omega * ((vy * cos(L) * cos(AZ)) - vx_w * sin(L))
    Mach = v/vs
    
    """ initializing variables for plots"""
    t_list = [t]
    x_list = [x]
    y_list = [y]
    vy_list = [vy]
    v_list = [v]
    phi0_list = [phi0]
    Mach_list = [Mach]
    
    while y >= 0:
        phi0 = phi0
        """drag calculation with variable density, Temp and sound speed"""
        rho = rho0 * ((Tb / (Tb + Lb * (y - y0))) ** (1 + (g * M / (R *Lb))))
        vs = 20.05 * np.sqrt(Tb + Lb * (y - y0)) # m/s speed of sound as a function of temperature
        Cd3 = K3 / sqrt(v / vs)
        Mach = v/vs
        
        """Determining drag regime"""
        if v > 1.2 * vs: #supersonic
            Cd = Cd3
        elif v < 0.8 * vs: #subsonic
            Cd = Cd1
        else: #transonic
            Cd = ((Cd3 - Cd1)*(v/vs - 0.8)/(0.4)) + Cd1
            
        """Acceleration due to Coriolis"""
        ax_c = -2*omega*((vz_w*sin(L))+ vy*cos(L)*sin(AZ))
        ay_c = 2*omega*((vz_w*cos(L)*cos(AZ))+ vx_w*cos(L)*sin(AZ))
        az_c = -2*omega*((vy*cos(L)*cos(AZ))- vx_w*sin(L))

        """Total acceleration calcs"""
        if vx > 0:
            ax = -0.5*rho*((vx-wx)**2)*Cd*Area/m + ax_c
        else:
            ax = 0

        """ Vy before and after peak"""
        if vy > 0:
            ay = (-0.5 * rho * (vy ** 2) * Cd * Area / m) - g + ay_c
        else:
            ay = (0.5 * rho * (vy ** 2) * Cd * Area / m) - g + ay_c
            az = az_c
            vx = vx + ax*dt # vx without wind
            # vx_w = vx with drag and no wind + wind
            vx_w = vx + 2*wx*(1-(vx/v0*np.cos(radians(phi0))))
            vy = vy + ay*dt
            vz = vz + az*dt
            vz_w = vz + wz*(1-(vx/v0*np.cos(radians(phi0))))
            """projectile velocity"""
            v = sqrt(vx_w**2 + vy**2 + vz**2)
            """new x, y, z positions"""
            x = x + vx_w*dt
            y = y + vy*dt
            z = z + vz_w*dt
            if y_max <= y:
                y_max = y
                phi_impact = degrees(atan(vy/vx)) #impact angle in degrees
                """ appends selected data for ability to plot"""
                t_list.append(t)
                x_list.append(x)
                y_list.append(y)
                vy_list.append(vy)
                v_list.append(v)
                phi0_list.append(phi0)
                Mach_list.append(Mach)
                if y < 0:
                    break
                    t += dt
                    
            """Range table output"""
            table[i] = ('%.f' % phi0, '%.3f' % phi_impact, '%.1f' % x,'%.2f' % z, '%.1f' % y_max, '%.1f' % vx_w,'%.1f' % vz,'%.1f' % vy,'%.1f' % v,'%.2f' %Mach, '%.1f' % t)

    """ Plot"""
    plt.plot(x_list, y_list, label='%d°' % phi0)#plt.plot(x_list, y_list, label='%d°' % phi0)
    plt.title('Altitude versus Range')
    plt.ylabel('Altitude (m)')
    plt.xlabel('Range (m)')
    plt.axis([0, 30000, 0, 15000])
    plt.grid(True)
    
print(table)

legend = plt.legend(title="Firing Angle",loc=0, fontsize='small', fancybox=True)

plt.show()

Thank you in advance

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

在梵高的星空下 2025-01-17 20:12:31

我应该使用哪个编辑器？

就我个人而言，我更喜欢VSCode，但是Sublime 也很受欢迎。如果你真的想使用准系统，请尝试 Vim。这三个都是完全免费的。

代码错误

扫描代码片段后，您似乎陷入了无限循环，您使用语句 while y >= 0 输入该循环。当您按下 Ctrl+C 时，您总是会看到第 103 行，原因很可能是因为这需要最长的时间，因此更有可能在任何给定时间到达该处。

请注意，目前，您只能通过此分支转义 while 循环：

    if y_max <= y:

        y_max= y
        phi_impact = degrees(atan(vy/vx)) #impact angle in degrees
        """ appends selected data for ability to plot"""
        t_list.append(t)
        x_list.append(x)
        y_list.append(y)
        vy_list.append(vy)
        v_list.append(v)
        phi0_list.append(phi0)
        Mach_list.append(Mach)
        if y < 0:
            break
            t += dt

这意味着如果 ymax 永远不会低于 y，或 y 永远不会低于零，那么你将无限循环。当然，我没有深入研究您的代码，但从表面上看，y_max 永远不会递减（这意味着它始终至少等于 y >）。此外，y 仅当您执行 y = y + vy*dt 时才会更新，这只会增加 y > if vy >= 0 （我假设 dt 始终为正）。

调试

正如 @Giacomo Catenazzi 建议的那样，尝试在 while 循环顶部打印 y 和 y_max 并查看它们在代码运行时如何变化。我怀疑它们并没有像你预期的那样减少。

Which Editor Should I Use?

Personally, I prefer VSCode, but Sublime is also pretty popular. If you really want to go barebones, try Vim. All three are completely free.

Code Errors

After scanning your code snippet, it appears that you are caught in an infinite loop, which you enter with the statement while y >= 0. The reason you always get line 103 when you hit Ctrl+C is likely because that takes the longest, making it more likely to land there at any given time.

Note that currently, you can only escape your while loop through this branch:

    if y_max <= y:

        y_max= y
        phi_impact = degrees(atan(vy/vx)) #impact angle in degrees
        """ appends selected data for ability to plot"""
        t_list.append(t)
        x_list.append(x)
        y_list.append(y)
        vy_list.append(vy)
        v_list.append(v)
        phi0_list.append(phi0)
        Mach_list.append(Mach)
        if y < 0:
            break
            t += dt

This means that if ymax never drops below y, or y never drops below zero, then you will infinitely loop. Granted, I haven't looked at your code in any great depth, but from the surface it appears that y_max is never decremented (meaning it will always be at least equal to y). Furthermore, y is only updated when you do y = y + vy*dt, which will only ever increase y if vy >= 0 (I assume dt is always positive).

Debugging

As @Giacomo Catenazzi suggested, try printing out y and y_max at the top of the while loop and see how they change as your code runs. I suspect they are not decrementing like you expected.

回复收藏 0 原文

~没有更多了~