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
