mirror of
https://gitlab.astro-wise.org/micado/atmosphericmodel.git
synced 2026-10-30 05:01:02 +00:00
233 lines
8.0 KiB
Python
233 lines
8.0 KiB
Python
#!/usr/bin/env python
|
|
|
|
""" Tables.py - Make tables of atmospheric properties (Python 3)
|
|
|
|
Adapted by
|
|
Richard J. Kwan, Lightsaber Computing
|
|
from original programs by
|
|
Ralph L. Carmichael, Public Domain Aeronautical Software
|
|
|
|
Revision History
|
|
Date Vers Person Statement of Changes
|
|
2004 Oct 04 1.0 RJK Initial program
|
|
2017 Jun 04 1.1 RLC All indents are with spaces; all prints in ( )
|
|
New version for Python 3 that does integer div with //
|
|
"""
|
|
|
|
import sys, math
|
|
|
|
version = "1.1 (2017 Jun 04)"
|
|
greeting = "Tables - A Python program to compute atmosphere tables"
|
|
author = "Ralph L. Carmichael, Public Domain Aeronautical Software"
|
|
modifier = "Richard J. Kwan, Lightsaber Computing"
|
|
farewell = "Four files added to your directory."
|
|
finalmess = "Normal termination of tables."
|
|
|
|
# P H Y S I C A L C O N S T A N T S
|
|
|
|
FT2METERS = 0.3048 # mult. ft. to get meters (exact)
|
|
KELVIN2RANKINE = 1.8 # mult deg K to get deg R
|
|
PSF2NSM = 47.880258 # mult lb/sq.ft to get sq.m
|
|
SCF2KCM = 515.379 # mult slugs/cu.ft to get kg/cu.m
|
|
TZERO = 288.15 # sea-level temperature, kelvins
|
|
PZERO = 101325.0 # sea-level pressure, N/sq.m
|
|
RHOZERO = 1.225 # sea-level density, kg/cu.m
|
|
AZERO = 340.294 # speed of sound at S.L. m/sec
|
|
BETAVISC = 1.458E-6 # viscosity constant
|
|
SUTH = 110.4 # Sutherland's constant, kelvins
|
|
|
|
def LongUSTable():
|
|
Itxt = open('us1py.prt', 'w')
|
|
Itxt.write(' alt sigma delta theta ')
|
|
Itxt.write(' temp press dens a visc k.visc\n')
|
|
Itxt.write(' Kft ')
|
|
Itxt.write(' degR lb/sq.ft s/cu.ft fps s/ft-s sq.ft/s\n')
|
|
|
|
for i in range(-1, 57):
|
|
altKm=5*i*FT2METERS
|
|
(sigma, delta, theta) = Atmosphere(altKm)
|
|
Itxt.write("%4i %9.3E %9.3E %6.4f " %
|
|
(5*i, sigma, delta, theta))
|
|
temp=(KELVIN2RANKINE*TZERO)*theta
|
|
pressure=(PZERO/PSF2NSM)*delta
|
|
density=(RHOZERO/SCF2KCM)*sigma
|
|
asound=(AZERO/FT2METERS)*math.sqrt(theta)
|
|
Itxt.write("%5.1f %8.3E %9.3E %6.1f" %
|
|
(temp, pressure, density, asound))
|
|
viscosity=(1.0/PSF2NSM)*MetricViscosity(theta)
|
|
kinematicViscosity=viscosity/density
|
|
Itxt.write("%6.3f %8.2E\n" %
|
|
(1.0E6*viscosity, kinematicViscosity))
|
|
|
|
Itxt.close()
|
|
|
|
def ShortUSTable():
|
|
Itxt = open('us2py.prt', 'w')
|
|
Itxt.write(' alt sigma delta theta ')
|
|
Itxt.write(' temp press dens a visc k.visc ratio\n')
|
|
Itxt.write(' Kft ')
|
|
Itxt.write(' degR psf s/cu.ft fps s/ft-sec sq.ft/s 1/ft\n')
|
|
|
|
for i in range(-1, 66):
|
|
altKm=i*FT2METERS
|
|
(sigma, delta, theta) = SimpleAtmosphere(altKm)
|
|
Itxt.write("%4i %6.4f %6.4f %6.4f" %
|
|
(i, sigma, delta, theta))
|
|
temp=KELVIN2RANKINE*TZERO*theta
|
|
pressure=PZERO*delta/47.88
|
|
density=RHOZERO*sigma/515.379
|
|
asound=(AZERO/FT2METERS)*math.sqrt(theta);
|
|
Itxt.write("%6.1f %6.1f %9.7f %6.1f" %
|
|
(temp, pressure, density, asound))
|
|
viscosity=(1.0/PSF2NSM)*MetricViscosity(theta)
|
|
kinematicViscosity=viscosity/density
|
|
vratio=asound/kinematicViscosity
|
|
Itxt.write("%6.3f %8.2E %4.2f\n"%
|
|
(1.0E6*viscosity, kinematicViscosity, 1.0E-6*vratio))
|
|
|
|
Itxt.close()
|
|
|
|
def LongSITable():
|
|
Itxt = open('si1py.prt', 'w')
|
|
Itxt.write(' alt sigma delta theta ')
|
|
Itxt.write(' temp press dens a visc k.visc\n')
|
|
Itxt.write(' Km ')
|
|
Itxt.write(' K N/sq.m kg/cu.m m/sec kg/m-s sq.m/s\n')
|
|
|
|
for i in range(-4,173):
|
|
altKm=0.5*i
|
|
(sigma, delta, theta) = Atmosphere(altKm)
|
|
Itxt.write("%4.1f %8.4E %8.4E %6.4f" %
|
|
(0.5*i, sigma,delta,theta))
|
|
temp=TZERO*theta
|
|
pressure=PZERO*delta
|
|
density=RHOZERO*sigma
|
|
asound=AZERO*math.sqrt(theta)
|
|
Itxt.write("%6.1f %8.5E %8.3E %5.1f" %
|
|
(temp,pressure,density,asound))
|
|
viscosity=MetricViscosity(theta)
|
|
kinematicViscosity=viscosity/density
|
|
Itxt.write("%6.2f %8.2E\n" %
|
|
(1.0E6*viscosity, kinematicViscosity))
|
|
|
|
Itxt.close()
|
|
|
|
def ShortSITable():
|
|
Itxt = open('si2py.prt', 'w')
|
|
Itxt.write(' alt sigma delta theta ')
|
|
Itxt.write(' temp press dens a visc k.visc ratio\n')
|
|
Itxt.write(' Km ')
|
|
Itxt.write(' K N/sq.m kcm m/sec kg/m-s sq.m/s 1/m\n')
|
|
|
|
for i in range(-1, 41):
|
|
altKm=0.5*i
|
|
(sigma, delta, theta) = SimpleAtmosphere(altKm)
|
|
Itxt.write("%4.1f %6.4f %6.4f %6.4f" %
|
|
(altKm, sigma, delta, theta))
|
|
temp=TZERO*theta
|
|
pressure=PZERO*delta
|
|
density=RHOZERO*sigma
|
|
asound=AZERO*math.sqrt(theta)
|
|
Itxt.write("%6.1f %6.0f %5.3f %5.1f" %
|
|
(temp, pressure,density,asound))
|
|
viscosity=MetricViscosity(theta)
|
|
kinematicViscosity=viscosity/density
|
|
vratio=asound/kinematicViscosity
|
|
Itxt.write("%6.2f %8.2E %5.2f\n" %
|
|
(1.0E6*viscosity, kinematicViscosity, 1.0E-6*vratio))
|
|
|
|
Itxt.close()
|
|
|
|
def Atmosphere(alt):
|
|
""" Compute temperature, density, and pressure in standard atmosphere.
|
|
Correct to 86 km. Only approximate thereafter.
|
|
Input:
|
|
alt geometric altitude, km.
|
|
Return: (sigma, delta, theta)
|
|
sigma density/sea-level standard density
|
|
delta pressure/sea-level standard pressure
|
|
theta temperature/sea-level std. temperature
|
|
"""
|
|
|
|
REARTH = 6369.0 # radius of the Earth (km)
|
|
GMR = 34.163195
|
|
NTAB = 8 # length of tables
|
|
|
|
htab = [ 0.0, 11.0, 20.0, 32.0, 47.0,
|
|
51.0, 71.0, 84.852 ]
|
|
ttab = [ 288.15, 216.65, 216.65, 228.65, 270.65,
|
|
270.65, 214.65, 186.946 ]
|
|
ptab = [ 1.0, 2.2336110E-1, 5.4032950E-2, 8.5666784E-3, 1.0945601E-3,
|
|
6.6063531E-4, 3.9046834E-5, 3.68501E-6 ]
|
|
gtab = [ -6.5, 0.0, 1.0, 2.8, 0, -2.8, -2.0, 0.0 ]
|
|
|
|
h = alt*REARTH/(alt+REARTH) # geometric to geopotential altitude
|
|
|
|
i=0; j=len(htab)-1
|
|
while (j > i+1):
|
|
k = (i+j)//2 # this is floor division in Python 3
|
|
if h < htab[k]:
|
|
j = k
|
|
else:
|
|
i = k
|
|
tgrad = gtab[i] # temp. gradient of local layer
|
|
tbase = ttab[i] # base temp. of local layer
|
|
deltah=h-htab[i] # height above local base
|
|
tlocal=tbase+tgrad*deltah # local temperature
|
|
theta = tlocal/ttab[0] # temperature ratio
|
|
|
|
if 0.0 == tgrad:
|
|
delta=ptab[i]*math.exp(-GMR*deltah/tbase)
|
|
else:
|
|
delta=ptab[i]*math.pow(tbase/tlocal, GMR/tgrad)
|
|
sigma = delta/theta
|
|
return ( sigma, delta, theta )
|
|
|
|
def SimpleAtmosphere(alt):
|
|
""" Compute temperature, density, and pressure in simplified
|
|
standard atmosphere.
|
|
|
|
Correct to 20 km. Only approximate thereafter.
|
|
|
|
Input:
|
|
alt geometric altitude, km.
|
|
Return: (sigma, delta, theta)
|
|
sigma density/sea-level standard density
|
|
delta pressure/sea-level standard pressure
|
|
theta temperature/sea-level std. temperature
|
|
"""
|
|
REARTH = 6369.0 # radius of the Earth (km)
|
|
GMR = 34.163195 # gas constant
|
|
|
|
h = alt*REARTH/(alt+REARTH) # geometric to geopotential altitude
|
|
|
|
if h<11.0: # troposphere
|
|
theta=(288.15-6.5*h)/288.15
|
|
delta=math.pow(theta, GMR/6.5)
|
|
else: # stratosphere
|
|
theta=216.65/288.15
|
|
delta=0.2233611*math.exp(-GMR*(h-11.0)/216.65)
|
|
sigma = delta/theta
|
|
return ( sigma, delta, theta )
|
|
|
|
def MetricViscosity(theta):
|
|
t=theta*TZERO
|
|
return BETAVISC*math.sqrt(t*t*t)/(t+SUTH)
|
|
|
|
def main():
|
|
print ("Executing ", sys.argv[0])
|
|
print (greeting)
|
|
print (author)
|
|
if modifier != "":
|
|
print ("Modified by ", modifier)
|
|
print (" version ", version)
|
|
LongUSTable()
|
|
ShortUSTable()
|
|
LongSITable()
|
|
ShortSITable()
|
|
|
|
print (farewell)
|
|
print (finalmess)
|
|
|
|
main()
|