Advertisement
meithan

US1976.py

Sep 5th, 2016
409
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 10.79 KB | None | 0 0
  1. # -*- coding:utf-8 -*-
  2. from math import exp, sqrt
  3.  
  4. ################################################################################
  5.  
  6. # This class implements the U.S. Standard Atmosphere 1976 (US1976) atmospheric
  7. # model up to 86 km of altitude.
  8.  
  9. # Under this model, the atmosphere below 86 km is approximated as 7 distinct
  10. # layers in which temperature is either constant or varies linearly with height,
  11. # having uniform chemical composition throughout.
  12.  
  13. # Above 86 km, the chemical composition is not uniform and depends on photo-
  14. # chemical processes and molecular diffusion. This implementation does not
  15. # model conditions above this altitude and will return density and pressure of
  16. # zero, and leave temperature and temperature-dependent quantities undefined.
  17. # For negative altitudes, the model returns sea-level values.
  18.  
  19. # The main API comprises the following class methods, which all take as input
  20. # the (geometric) altitude, in meters, and return the corresponding atmospheric
  21. # variable, in SI units:
  22. #  * density(altitude)
  23. #  * pressure(altitude)
  24. #  * temperature(altitude)
  25. #  * sound_speed(altitude)
  26. #  * viscosity(altitude)
  27.  
  28. # Alternatively, one can directly call the values(altitude, var=None) method,
  29. # passing either no var argument or var="all", which will return a tuple with
  30. # either (density, pressure, temperature) or (density, pressure, temperature,
  31. # sound speed, viscosity), respectively. This is faster than calling the above
  32. # methods individually when one wants multiple variables at the same height.
  33.  
  34. # Reference document for the U.S. Standard Atmosphere 1976:
  35. # http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19770009539.pdf
  36.  
  37. # A summary of the definition can be found here:
  38. # http://www.pdas.com/atmosdef.html
  39.  
  40. # Copyright (C) 2016 Meithan West (aka Juan C. Toledo)
  41.  
  42. # This program is free software: you can redistribute it and/or modify
  43. # it under the terms of the GNU General Public License as published by
  44. # the Free Software Foundation, either version 3 of the License, or
  45. # (at your option) any later version.
  46.  
  47. # This program is distributed in the hope that it will be useful,
  48. # but WITHOUT ANY WARRANTY; without even the implied warranty of
  49. # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  50. # GNU General Public License for more details.
  51.  
  52. # You should have received a copy of the GNU General Public License
  53. # along with this program. If not, see http://www.gnu.org/licenses/.
  54.  
  55. ################################################################################
  56.  
  57. class US1976:
  58.  
  59.   # Initialize model
  60.   def __init__(self):
  61.  
  62.     # Layers of the model
  63.     # Each layer is specified by a tuple of the form
  64.     #   (base_alt, base_pres, base_temp, lapse_rate)
  65.     # where
  66.     #   base_alt: *geometric* altitude at base of layer, in m
  67.     #   base_pres: pressure at base of layer, in Pa
  68.     #   base_temp: temperature at base of layer, in K
  69.     #   lapse_rate: temperature lapse rate in the layer, in K/km'
  70.     #   NOTE: temperature lapse rates are specified in Kelvins per km of
  71.     #   *geopotential* altitude (not geometric altitude).
  72.     # Each layer extends up to the base altitude of the next (with the
  73.     # last extending up to zmax, defined below).
  74.     # The base altitudes correspond to 0, 11, 20, 32, 47, 51 and 71 km
  75.     # of geopotential altitude, respectively.
  76.     self.layers = (
  77.     (0.0,     101325.0, 288.15, -6.5),
  78.     (11019.0, 22632.0,  216.65,  0.0),
  79.     (20063.1, 5474.9,   216.65, +1.0),
  80.     (32161.9, 868.02,   228.65, +2.8),
  81.     (47350.1, 110.91,   270.65,  0.0),
  82.     (51412.5, 66.939,   270.65, -2.8),
  83.     (71802.0, 3.9564,   214.65, -2.0))
  84.     self.nlayers = len(self.layers)
  85.  
  86.     # Maximum (geometric) altitude of the model, m
  87.     self.zmax = 86000.0
  88.  
  89.     # Universal gas constant, J/(K mol)
  90.     self.Rg = 8.31432
  91.  
  92.     # Molecular weight of dry air, kg/mol
  93.     self.M_air = 28.9644e-3
  94.  
  95.     # Ratio of specific heats of air
  96.     self.gamma = 1.4
  97.  
  98.     # Beta constant for viscosity, kg/(s m K^1/2)
  99.     self.beta = 1.458e-6
  100.  
  101.     # Sutherland's constant, K
  102.     self.S = 110.4
  103.  
  104.     # Radius of the Earth (average for 45°N latitude), m
  105.     self.RE = 6356.766e3
  106.  
  107.     # Standard gravity, m/s^2
  108.     self.g0 = 9.80665
  109.  
  110.     # Specific gas constant for dry air, J/(K kg)
  111.     self.R_air = self.Rg/self.M_air
  112.  
  113.     # Sea-level (zero altitude) values, SI units
  114.     self.P_sl = self.layers[0][1]
  115.     self.T_sl = self.layers[0][2]
  116.     self.rho_sl = self.P_sl/(self.R_air*self.T_sl)
  117.     self.cs_sl = sqrt(self.gamma * self.R_air * self.T_sl)
  118.     self.mu_sl = self.beta*self.T_sl**(3./2.)/(self.T_sl+self.S)
  119.  
  120.   # ==================================
  121.  
  122.   # Returns the index of the layer the given altitude lies in
  123.   # (or None upon error); altitude is geometric and in meters
  124.   def get_layer(self, altitude):
  125.  
  126.     if altitude < 0.0 or altitude > self.zmax:
  127.       return None
  128.  
  129.     # Binary search for layer
  130.     low = 0
  131.     high = len(self.layers)-1
  132.     while low != high:
  133.       mid = (low+high)//2
  134.       if altitude > self.layers[mid+1][0]:
  135.         low = mid+1
  136.       else:
  137.         high = mid
  138.     return low
  139.  
  140.   # ==================================
  141.  
  142.   # Converts geometric altitude to geopotential altitude
  143.   def to_geopot_alt(self, geom_alt):
  144.     return geom_alt*self.RE / (self.RE + geom_alt)
  145.  
  146.   # Converts geopotential altitude to geometric altitude
  147.   def to_geomet_alt(self, geop_alt):
  148.     return geop_alt*self.RE / (self.RE - geop_alt)
  149.  
  150.   # ==================================
  151.  
  152.   # Definitions of sound speed and viscosity as functions of T
  153.  
  154.   def sound_speed_from_T(self, T):
  155.     return sqrt(self.gamma * self.R_air * T)
  156.  
  157.   def viscosity_from_T(self, T):    
  158.     return self.beta*T**(3./2.)/(T+self.S)
  159.  
  160.   # ==================================
  161.  
  162.   # API functions
  163.  
  164.   # The following functions are wrappers that compute and return a single
  165.   # atmospheric variable when given a geometric altitude
  166.  
  167.   def density(self, altitude):
  168.     """Compute atmospheric density at altitude.  
  169.    Input:
  170.      altitude: geometric altitude, in m.
  171.    Output:
  172.      A single float with density, in kg/m^3.
  173.    """
  174.     return self.values(altitude, var="dens")
  175.  
  176.   def pressure(self, altitude):
  177.     """Compute atmospheric pressure at altitude.  
  178.    Input:
  179.      altitude: geometric altitude, in m.
  180.    Output:
  181.      A single float with pressure, in Pa.
  182.    """
  183.     return self.values(altitude, var="pres")
  184.  
  185.   def temperature(self, altitude):
  186.     """Compute atmospheric temperature at altitude.  
  187.    Input:
  188.      altitude: geometric altitude, in m.
  189.    Output:
  190.      A single float with temperature, in K.
  191.    """  
  192.     return self.values(altitude, var="temp")
  193.  
  194.   def sound_speed(self, altitude):
  195.     """Compute atmospheric sound speed at altitude.  
  196.    Input:
  197.      altitude: geometric altitude, in m.
  198.    Output:
  199.      A single float with sound speed, in m/s:
  200.    """
  201.     return self.sound_speed_from_T(self.temperature(altitude))
  202.    
  203.   def viscosity(self, altitude):
  204.     """Compute atmospheric dynamic viscosity at altitude.  
  205.    Input:
  206.      altitude: geometric altitude, in m.
  207.    Output:
  208.      A single float with dynamic viscosity, Pa s.
  209.    """
  210.     return self.viscosity_from_T(self.temperature(altitude))
  211.  
  212.   # The following function does the actual computation of atmospheric values
  213.   def values(self, altitude, var=None):
  214.     """Compute atmospheric values at altitude
  215.  
  216.    Input:
  217.     altitude: geometric altitude, in m.
  218.  
  219.    Output:
  220.      Depends on the value of argument 'var':
  221.      None: a (density, pressure, temperature) tuple; this is the default
  222.      "all": a (density, pressure, temperature, sound speed, viscosity) tuple
  223.      "dens": a float with density, in kg/m^3
  224.      "pres": a float with pressure, in Pa
  225.      "temp": a float with temperature, in K
  226.      "sound": a float with the speed of sound, in m/s
  227.      "visc": a float with the dynamic viscosity, in Pa s
  228.  
  229.    If the provided altitude is above zmax=86 km it will return zero values for
  230.    density and pressure and an undefined temperature. For negative altitudes,
  231.    it will return the sea-level values.
  232.    """
  233.     if altitude > self.zmax:
  234.  
  235.       if var == 'pres' or var == 'dens':
  236.         return 0.0
  237.       elif var == 'temp' or var == "sound" or var == "visc":
  238.         return float('NaN')
  239.       elif var == None:
  240.         return (0.0, 0.0, float('NaN'))
  241.       else:
  242.         return (0.0, 0.0, float('NaN'), float('NaN'), float('NaN'))
  243.        
  244.     elif altitude <= 0:
  245.  
  246.       if var == 'dens': return self.rho_sl
  247.       elif var == 'pres': return self.P_sl
  248.       elif var == 'temp': return self.T_sl
  249.       elif var == 'sound': return self.cs_sl
  250.       elif var == 'visc': return self.mu_sl
  251.       elif var == None: return (self.rho_sl, self.P_sl, self.T_sl)
  252.       elif var == 'all': return (self.rho_sl, self.P_sl, self.T_sl,\
  253.       self.cs_sl, self.mu_sl)
  254.  
  255.     else:
  256.  
  257.       # Determine in which layer the given altitude is
  258.       layer = self.get_layer(altitude)
  259.  
  260.       # Values at layer base and layer lapse rate
  261.       Z0 = self.layers[layer][0]
  262.       P0 = self.layers[layer][1]
  263.       T0 = self.layers[layer][2]
  264.       TL = self.layers[layer][3]
  265.       TL = TL/1.0e3    # per meter
  266.      
  267.       # Convert geometric altitudes to geopotential ones
  268.       H0 = self.to_geopot_alt(Z0)
  269.       H = self.to_geopot_alt(altitude)
  270.      
  271.       # Compute temperature
  272.       if TL == 0:
  273.         temp = T0
  274.       else:
  275.         temp = T0 + (H-H0)*TL
  276.       if var == "temp": return temp
  277.      
  278.       # Compute speed of sound
  279.       if var == "all" or var == "sound":
  280.         sound = self.sound_speed_from_T(temp)
  281.         if var == "sound": return sound
  282.        
  283.       # Compute dynamic viscosity
  284.       if var == "all" or var == "visc":
  285.         visc = self.viscosity_from_T(temp)
  286.         if var == "visc": return visc
  287.      
  288.       # Compute pressure
  289.       if TL == 0: pres = P0 * exp(-(H-H0)*self.g0/(self.R_air*T0))
  290.       else: pres = P0 * pow(temp/T0,-self.g0/(self.R_air*TL))
  291.       if var == "pres": return pres
  292.      
  293.       # Compute density
  294.       dens = pres/(self.R_air*temp)
  295.       if var == "dens": return dens
  296.  
  297.       # Return tuples with multiple values
  298.       if var == None:
  299.         return (dens, pres, temp)
  300.       else:
  301.         return (dens, pres, temp, sound, visc)
  302.  
  303. ################################################################################
  304.  
  305. # Unit testing: print a table with values every 1 km
  306. if __name__ == "__main__":
  307.   atmo = US1976()
  308.   z = 0.0
  309.   print("Altitude (geometric), altitude (geopotential), density, pressure, \
  310. temperature, sound speed, viscosity, layer #")
  311.   while (z <= 86000.0):
  312.     dens, pres, temp, cs, mu = atmo.values(z,var='all')
  313.     print("%6.3f, %6.3f, %.8f, %10.3f, %.2f, %.1f, %.2f, %i" %
  314.     (z/1e3, atmo.to_geopot_alt(z)/1e3, dens, pres, temp, cs, mu*1e6, atmo.get_layer(z)+1))
  315.     z += 1000.0
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement