US1976.py

# -*- coding:utf-8 -*-
from math import exp, sqrt

################################################################################

# This class implements the U.S. Standard Atmosphere 1976 (US1976) atmospheric
# model up to 86 km of altitude.

# Under this model, the atmosphere below 86 km is approximated as 7 distinct
# layers in which temperature is either constant or varies linearly with height,
# having uniform chemical composition throughout.

# Above 86 km, the chemical composition is not uniform and depends on photo-
# chemical processes and molecular diffusion. This implementation does not
# model conditions above this altitude and will return density and pressure of
# zero, and leave temperature and temperature-dependent quantities undefined.
# For negative altitudes, the model returns sea-level values.

# The main API comprises the following class methods, which all take as input
# the (geometric) altitude, in meters, and return the corresponding atmospheric
# variable, in SI units:
#  * density(altitude)
#  * pressure(altitude)
#  * temperature(altitude)
#  * sound_speed(altitude)
#  * viscosity(altitude)

# Alternatively, one can directly call the values(altitude, var=None) method,
# passing either no var argument or var="all", which will return a tuple with
# either (density, pressure, temperature) or (density, pressure, temperature,
# sound speed, viscosity), respectively. This is faster than calling the above
# methods individually when one wants multiple variables at the same height.

# Reference document for the U.S. Standard Atmosphere 1976:
# http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19770009539.pdf

# A summary of the definition can be found here:
# http://www.pdas.com/atmosdef.html

# Copyright (C) 2016 Meithan West (aka Juan C. Toledo)

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program. If not, see http://www.gnu.org/licenses/.

################################################################################

class US1976:

  # Initialize model
  def __init__(self):

    # Layers of the model
    # Each layer is specified by a tuple of the form
    #   (base_alt, base_pres, base_temp, lapse_rate)
    # where
    #   base_alt: *geometric* altitude at base of layer, in m
    #   base_pres: pressure at base of layer, in Pa
    #   base_temp: temperature at base of layer, in K
    #   lapse_rate: temperature lapse rate in the layer, in K/km'
    #   NOTE: temperature lapse rates are specified in Kelvins per km of
    #   *geopotential* altitude (not geometric altitude).
    # Each layer extends up to the base altitude of the next (with the
    # last extending up to zmax, defined below).
    # The base altitudes correspond to 0, 11, 20, 32, 47, 51 and 71 km
    # of geopotential altitude, respectively.
    self.layers = (
    (0.0,     101325.0, 288.15, -6.5),
    (11019.0, 22632.0,  216.65,  0.0),
    (20063.1, 5474.9,   216.65, +1.0),
    (32161.9, 868.02,   228.65, +2.8),
    (47350.1, 110.91,   270.65,  0.0),
    (51412.5, 66.939,   270.65, -2.8),
    (71802.0, 3.9564,   214.65, -2.0))
    self.nlayers = len(self.layers)

    # Maximum (geometric) altitude of the model, m
    self.zmax = 86000.0

    # Universal gas constant, J/(K mol)
    self.Rg = 8.31432

    # Molecular weight of dry air, kg/mol
    self.M_air = 28.9644e-3

    # Ratio of specific heats of air
    self.gamma = 1.4

    # Beta constant for viscosity, kg/(s m K^1/2)
    self.beta = 1.458e-6

    # Sutherland's constant, K
    self.S = 110.4

    # Radius of the Earth (average for 45°N latitude), m
    self.RE = 6356.766e3

    # Standard gravity, m/s^2
    self.g0 = 9.80665

    # Specific gas constant for dry air, J/(K kg)
    self.R_air = self.Rg/self.M_air

    # Sea-level (zero altitude) values, SI units
    self.P_sl = self.layers[0][1]
    self.T_sl = self.layers[0][2]
    self.rho_sl = self.P_sl/(self.R_air*self.T_sl)
    self.cs_sl = sqrt(self.gamma * self.R_air * self.T_sl)
    self.mu_sl = self.beta*self.T_sl**(3./2.)/(self.T_sl+self.S)

  # ==================================

  # Returns the index of the layer the given altitude lies in
  # (or None upon error); altitude is geometric and in meters
  def get_layer(self, altitude):

    if altitude < 0.0 or altitude > self.zmax:
      return None

    # Binary search for layer
    low = 0
    high = len(self.layers)-1
    while low != high:
      mid = (low+high)//2
      if altitude > self.layers[mid+1][0]:
        low = mid+1
      else:
        high = mid
    return low

  # ==================================

  # Converts geometric altitude to geopotential altitude
  def to_geopot_alt(self, geom_alt):
    return geom_alt*self.RE / (self.RE + geom_alt)

  # Converts geopotential altitude to geometric altitude
  def to_geomet_alt(self, geop_alt):
    return geop_alt*self.RE / (self.RE - geop_alt)

  # ==================================

  # Definitions of sound speed and viscosity as functions of T

  def sound_speed_from_T(self, T):
    return sqrt(self.gamma * self.R_air * T)

  def viscosity_from_T(self, T):
    return self.beta*T**(3./2.)/(T+self.S)

  # ==================================

  # API functions

  # The following functions are wrappers that compute and return a single
  # atmospheric variable when given a geometric altitude

  def density(self, altitude):
    """Compute atmospheric density at altitude.
    Input:
      altitude: geometric altitude, in m.
    Output:
      A single float with density, in kg/m^3.
    """
    return self.values(altitude, var="dens")

  def pressure(self, altitude):
    """Compute atmospheric pressure at altitude.
    Input:
      altitude: geometric altitude, in m.
    Output:
      A single float with pressure, in Pa.
    """
    return self.values(altitude, var="pres")

  def temperature(self, altitude):
    """Compute atmospheric temperature at altitude.
    Input:
      altitude: geometric altitude, in m.
    Output:
      A single float with temperature, in K.
    """
    return self.values(altitude, var="temp")

  def sound_speed(self, altitude):
    """Compute atmospheric sound speed at altitude.
    Input:
      altitude: geometric altitude, in m.
    Output:
      A single float with sound speed, in m/s:
    """
    return self.sound_speed_from_T(self.temperature(altitude))

  def viscosity(self, altitude):
    """Compute atmospheric dynamic viscosity at altitude.
    Input:
      altitude: geometric altitude, in m.
    Output:
      A single float with dynamic viscosity, Pa s.
    """
    return self.viscosity_from_T(self.temperature(altitude))

  # The following function does the actual computation of atmospheric values
  def values(self, altitude, var=None):
    """Compute atmospheric values at altitude

    Input:
     altitude: geometric altitude, in m.

    Output:
      Depends on the value of argument 'var':
      None: a (density, pressure, temperature) tuple; this is the default
      "all": a (density, pressure, temperature, sound speed, viscosity) tuple
      "dens": a float with density, in kg/m^3
      "pres": a float with pressure, in Pa
      "temp": a float with temperature, in K
      "sound": a float with the speed of sound, in m/s
      "visc": a float with the dynamic viscosity, in Pa s

    If the provided altitude is above zmax=86 km it will return zero values for
    density and pressure and an undefined temperature. For negative altitudes,
    it will return the sea-level values.
    """
    if altitude > self.zmax:

      if var == 'pres' or var == 'dens':
        return 0.0
      elif var == 'temp' or var == "sound" or var == "visc":
        return float('NaN')
      elif var == None:
        return (0.0, 0.0, float('NaN'))
      else:
        return (0.0, 0.0, float('NaN'), float('NaN'), float('NaN'))

    elif altitude <= 0:

      if var == 'dens': return self.rho_sl
      elif var == 'pres': return self.P_sl
      elif var == 'temp': return self.T_sl
      elif var == 'sound': return self.cs_sl
      elif var == 'visc': return self.mu_sl
      elif var == None: return (self.rho_sl, self.P_sl, self.T_sl)
      elif var == 'all': return (self.rho_sl, self.P_sl, self.T_sl,\
      self.cs_sl, self.mu_sl)

    else:

      # Determine in which layer the given altitude is
      layer = self.get_layer(altitude)

      # Values at layer base and layer lapse rate
      Z0 = self.layers[layer][0]
      P0 = self.layers[layer][1]
      T0 = self.layers[layer][2]
      TL = self.layers[layer][3]
      TL = TL/1.0e3    # per meter

      # Convert geometric altitudes to geopotential ones
      H0 = self.to_geopot_alt(Z0)
      H = self.to_geopot_alt(altitude)

      # Compute temperature
      if TL == 0:
        temp = T0
      else:
        temp = T0 + (H-H0)*TL
      if var == "temp": return temp

      # Compute speed of sound
      if var == "all" or var == "sound":
        sound = self.sound_speed_from_T(temp)
        if var == "sound": return sound

      # Compute dynamic viscosity
      if var == "all" or var == "visc":
        visc = self.viscosity_from_T(temp)
        if var == "visc": return visc

      # Compute pressure
      if TL == 0: pres = P0 * exp(-(H-H0)*self.g0/(self.R_air*T0))
      else: pres = P0 * pow(temp/T0,-self.g0/(self.R_air*TL))
      if var == "pres": return pres

      # Compute density
      dens = pres/(self.R_air*temp)
      if var == "dens": return dens

      # Return tuples with multiple values
      if var == None:
        return (dens, pres, temp)
      else:
        return (dens, pres, temp, sound, visc)

################################################################################

# Unit testing: print a table with values every 1 km
if __name__ == "__main__":
  atmo = US1976()
  z = 0.0
  print("Altitude (geometric), altitude (geopotential), density, pressure, \
temperature, sound speed, viscosity, layer #")
  while (z <= 86000.0):
    dens, pres, temp, cs, mu = atmo.values(z,var='all')
    print("%6.3f, %6.3f, %.8f, %10.3f, %.2f, %.1f, %.2f, %i" %
    (z/1e3, atmo.to_geopot_alt(z)/1e3, dens, pres, temp, cs, mu*1e6, atmo.get_layer(z)+1))
    z += 1000.0