United States Standard Atmosphere

ussa.py

__all__ = ['usstdatmdf', 'usstdatmds', 'M0', 'R', 'Rs', 'g', '__version__']
__version__ = '0.1.0'

__doc__ = """
This GIST implements the US Standard Atmosphere for Temperature, Pressure,
and Density. The constants and equations are implemented according to the U.S.
Standard Atmosphere (1976) and are only implemented through 84km. The results
can be presented as a Pandas DataFrame or an xarray Dataset.

Requirements
============

- numpy
- Pandas
- Optional, xarray

Examples
========

The code below assumes you have downloaded ussa.py. It imports the library,
creates standard atmosphere objects, and plots the results.

    import ussa
    import numpy as np

    # Create a default standard atmosphere dataframe at 100m resolution
    stdatmdf = ussa.usstdatmdf()
    # Plot pressure
    stdatmdf['P'].plot()

    # Create a default standard atmosphere dataset with 10 rows and cols
    stdatmds = ussa.usstdatmds(ROW=np.arange(10), COL=np.arange(10))
    # Plot temperature
    stdatmds['T'].sel(ROW=0, COL=0).plot(y='z')


Inventory
=========

usaa.stdatmdf : function
  Returns a dataframe of the standard atmosphere

usaa.stdatmdf : function
  Returns a dataset of the standard atmosphere

usaa.stdatmbdf : function
  Returns base-level constants described by the USSA 1976 include altitude,
  lapse rate, temperature (calculated), and pressure (calculated)

ussa.g : float
  Gravitational acceleration

ussa.M0 : float
  Atmospheric molecular weight (weighted mole fraction weight)

ussa.R : float
  Ideal gas constant

ussa.Rs : float
  Specific ideal gas constant (R / M0)


References
==========
U.S. Standard Atmosphere, 1976, U.S. Government Printing Office, Washington, D.C., 1976

"""
import numpy as np
import pandas as pd


# Downloaded from Wikipedia 2023-07-13
# https://en.wikipedia.org/wiki/U.S._Standard_Atmosphere
# Confirmed from:
# U.S. Standard Atmosphere, 1976, U.S. Government Printing Office, Washington, D.C., 1976
# zb : Geopotential height in meters above mean sea level
# Tb : Standard temperature in Kelvin
# Pb : Base pressure
stdatmbdf = pd.DataFrame(dict(
  zb=np.array([0, 11, 20, 32, 47, 51, 71, 84.852]) * 1e3,
  Tb=np.array([288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65, 187.15]),
  # Pb=np.array([
  #   101325, 22632.1, 5474.89, 868.019, 110.906, 66.9389, 3.95642, np.nan
  # ])
  # Recalculated for precision
  Pb=np.array([
    1.01325000e+05, 2.26320640e+04, 5.47489738e+03, 8.68018896e+02,
    1.10906346e+02, 6.69386887e+01, 3.95642202e+00, 3.73853320e-01
  ])
), index=np.arange(8))
# Add Lapse rate (K/m)
stdatmbdf.loc[0:6, 'Lb'] = np.diff(stdatmbdf['Tb']) / np.diff(stdatmbdf['zb'])
# Define  useful constants
R = 8.31432
M0 = 0.0289644
Rs = R / M0
g = 9.80665


def usstdatmdf(zs=None, debug=False):
  """
  Arguments
  ---------
  zs : array-like
    Geopotential altitude levels to return. If None, 100m spacing will be used.
  debug : bool
    If True, include intermediate calculations

  Returns
  -------
  zdf : pandas.DataFrame
    Standard atmosphere P, T, and density at zs levels
  """
  if zs is None:
    zs = np.linspace(0, 84852 + 100, 100)
  else:
    zs = np.asarray(zs)
  # In case zs is a dataframe series
  stdbdf = stdatmbdf

  idx = np.maximum(0, (zs[:, None] > stdbdf['zb'].values[None, :]).sum(1) - 1)
  zdf = stdbdf.loc[idx]
  zdf['z'] = zs
  zdf['T'] = zdf.eval('Tb + (z - zb) * Lb')
  # U.S. Standard Atmosphere, 1976, U.S. Government Printing Office, Washington, D.C., 1976
  # Equations 33a and 33B
  zdf['Pfactor_fT'] = zdf.eval('(Tb / T)')**(g / Rs / zdf['Lb'])
  zdf['Pfactor_iT'] = np.exp(-g * zdf.eval('(z - zb) / Tb') / Rs)
  # Select 33b when lapse rate is 0 and 33a otherwise
  zdf['Pfactor'] = np.where(
      zdf['Lb'] == 0, zdf['Pfactor_iT'], zdf['Pfactor_fT']
  )
  zdf['P'] = zdf.eval('Pb * Pfactor')
  zdf['dens'] = zdf['P'] / Rs / zdf['T']
  zdf = zdf.set_index('z')
  if not debug:
    zdf = zdf[['T', 'P', 'dens']]
  return zdf

def usstdatmds(zs=None, **dims):
  """
  Arguments
  ---------
  zs : array-like
    Geopotential altitude levels to return. If None, 100m spacing will be used.
  dims : mappable
    Dimensions to be added.

  Returns
  -------
  zdf : xarray.Dataset
    Standard atmosphere P, T, and density at zs levels
  """
  oned = {k: 1 for k in dims}
  temp = usstdatmdf(zs).to_xarray().expand_dims(**oned)
  temp['P'].attrs.update(
      long_name='Pressure', units='Pascals'
  )
  temp['T'].attrs.update(
      long_name='Temperature', units='K'
  )
  temp['dens'].attrs.update(
      long_name='Density', units='kg/m3'
  )
  for k in dims:
    temp.coords[k] = np.array([1])
  newdims = ['z'] + list(dims)
  temp = temp.sel(**dims, method='nearest').transpose(*newdims)
  for k in dims:
    temp.coords[k] = dims[k]
  return temp

Here is a basic method of downloading and applying.

The code below can be copied into a terminal and run:

wget -N https://gist.githubusercontent.com/barronh/00359ed2c07c7759cba6cbff0d6ce0dd/raw/ussa.py

python -ic "
import ussa
import matplotlib.pyplot as plt


stddf = ussa.usstdatmdf()
z = stddf.index.get_level_values('z')

fig, axx = plt.subplots(1, 3, sharey=True)
axx[0].plot(stddf['T'], z)
axx[0].set(xlim=(100, 300), xlabel='Temperature [K]')
axx[1].plot(stddf['dens'], z)
axx[1].set(xlim=(0, 2), xlabel='Density [kg/m3]')
axx[2].plot(stddf['P'] / 1000, z)
axx[2].set(xlim=(0, 100), xlabel='Pressure [kPa]')
fig.savefig('ussa1976_P.png')
"

To run in a jupyter notebook, run two cells as described below.

!wget -N https://gist.githubusercontent.com/barronh/00359ed2c07c7759cba6cbff0d6ce0dd/raw/ussa.py

import ussa
import matplotlib.pyplot as plt


stddf = ussa.usstdatmdf()
z = stddf.index.get_level_values('z')

fig, axx = plt.subplots(1, 3, sharey=True)
axx[0].plot(stddf['T'], z)
axx[0].set(xlim=(100, 300), xlabel='Temperature [K]')
axx[1].plot(stddf['dens'], z)
axx[1].set(xlim=(0, 2), xlabel='Density [kg/m3]')
axx[2].plot(stddf['P'] / 1000, z)
axx[2].set(xlim=(0, 100), xlabel='Pressure [kPa]')
fig.savefig('ussa1976_P.png')