Where is my code hanging (in an infinite loop)?

Question

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

Answer 1

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.