89 lines
3.0 KiB
Python
89 lines
3.0 KiB
Python
'''
|
|
Because of the interest of aerospace engineers and atmospheric scientists in
|
|
conditions at much higher altitudes than those of the U.S. Standard Atmosphere
|
|
1976(USSA 1976), members of the U.S. Committee on Extension to the Standard
|
|
Atmosphere(COESA 1976) agreed to extend it to 1000 km.
|
|
'''
|
|
|
|
import numpy as np
|
|
from pkg_resources import resource_filename
|
|
|
|
from .ussa76 import ussa76
|
|
from ..utils import Const
|
|
from ..utils.utils import alt_conver,check_altitude
|
|
|
|
from ..class_atmos import ATMOS
|
|
|
|
def coesa76(alts, alt_type='geometric'):
|
|
'''
|
|
Implements the U.S. Committee on Extension to the Standard Atmosphere(COESA 1976).
|
|
|
|
Usage:
|
|
[rhos, Ts, Ps] = coesa76(h)
|
|
|
|
Inputs:
|
|
alts -> [float list/array] geometric or geopotentail altitudes, [km]
|
|
|
|
Outputs:
|
|
rhos -> [float array] densities at a set of given altitudes, [kg/m^3]
|
|
Ts -> [float] temperatures ..., [K]
|
|
Ps -> [float] pressures ..., [Pa]
|
|
|
|
Note: the geometric altitudes should be in [-0.611,1000] km, otherwise the output will be extrapolated for those input altitudes.
|
|
|
|
Reference:
|
|
U.S. Standard Atmosphere, 1976, U.S. Government Printing Office, Washington, D.C.
|
|
http://www.braeunig.us/space/atmmodel.htm#USSA1976
|
|
https://docs.poliastro.space/en/stable/index.html
|
|
'''
|
|
|
|
# Base altitude for the COESA 1976, [km].
|
|
zb = np.array([86, 91, 100, 110, 120, 150, 200, 300, 500, 750, np.inf])
|
|
|
|
# load the coefficients used to approximate density and pressure above 86km
|
|
data_path = resource_filename('pyatmos', 'data/')
|
|
data = np.load(data_path+'coesa76_coeffs.npz')
|
|
rho_coeffs,p_coeffs = data['rho'],data['p']
|
|
|
|
R0 = Const.R0 # volumetric radius for the Earth, [km]
|
|
|
|
# Get geometric and geopotential altitudes
|
|
zs,hs = alt_conver(alts, alt_type)
|
|
|
|
# Test if altitudes are inside valid range
|
|
check_altitude(zs,(-0.611,1e3),'warning')
|
|
|
|
inds = np.zeros_like(zs,dtype=int)
|
|
rhos,Ts,Ps = np.zeros((3,len(zs)))
|
|
|
|
j = 0
|
|
for z,h in zip(zs,hs):
|
|
|
|
if z <= zb[0]:
|
|
rho,T,P,Cs,eta,Kc = ussa76(h)
|
|
else:
|
|
if z > zb[0] and z <= zb[1]:
|
|
T = 186.8673
|
|
elif z > zb[1] and z <= zb[3]:
|
|
T = 263.1905 - 76.3232 * np.sqrt(1 - ((z - 91) / 19.9429) ** 2)
|
|
elif z > zb[3] and z <= zb[4]:
|
|
T = 240 + 12 * (z - 110)
|
|
else:
|
|
epsilon = (z - 120) * (R0 + 120) / (R0 + z)
|
|
T = 1e3 - 640 * np.exp(-0.01875 * epsilon)
|
|
|
|
ind = np.where((z - zb) >= 0)[0][-1]
|
|
|
|
# A 4th order polynomial is used to approximate density and pressure.
|
|
# This is directly taken from: http://www.braeunig.us/space/atmmodel.htm
|
|
poly_rho = np.poly1d(rho_coeffs[ind])
|
|
poly_p = np.poly1d(p_coeffs[ind])
|
|
rho = np.exp(poly_rho(z))
|
|
P = np.exp(poly_p(z))
|
|
|
|
rhos[j],Ts[j],Ps[j] = rho,T,P
|
|
j += 1
|
|
|
|
info = {'rho':rhos,'T':Ts,'P':Ps}
|
|
|
|
return ATMOS(info) |