#!/usr/bin/python -tt
#=======================================================================
#                        General Documentation

"""Single-function module.

   See function docstring for description.
"""

#-----------------------------------------------------------------------
#                       Additional Documentation
#
# RCS Revision Code:
#   $Id: rho.py,v 1.1 2004/05/01 00:01:56 jlin Exp $
#
# Modification History:
# - 30 Apr 2003:  Original by Johnny Lin, Computation Institute,
#   University of Chicago.  Passed passably reasonable tests.
#
# Notes:
# - Written for Python 2.2.
# - Module docstrings can be tested using the doctest module.  To
#   test, execute "python rho.py".
# - See import statements throughout for non-"built-in" packages 
#   and modules required.
#
# Copyright (c) 2004 by Johnny Lin.  For licensing, distribution 
# conditions, contact information, and additional documentation see
# the URL http://www.johnny-lin.com/py_pkgs/atmqty/doc/.
#=======================================================================




#---------------- Module General Import and Declarations ---------------

#- Set module version number to package version number, etc.:

import package_version
__version__ = package_version.version
__author__  = package_version.author
__date__    = package_version.date
__credits__ = package_version.credits
del package_version




#--------------------------- Public Function ---------------------------

def rho(p_in, Tv_in, missing=1e+20):
    """Calculate air density (including moisture effects).

    Method Arguments:
    * p_in:  Pressure [hPa].  Required.
    * Tv_in:  Virtual temperature [K].  Required.
    * missing:  If p_in and/or Tv_in has missing values, this is the 
      missing value value.  Floating point scalar.  Default is 1e+20.

    p_in and Tv_in are Numeric floating point arrays of any number of
    dimensions and size, as long as they are of the same size and
    shape (or one of them is of size 1).

    Output:
    * Air density [kg/m**3].  Numeric array of same dimensions and 
      size as input.  If there are any missing values in output, those 
      values are set to the value in argument missing from the input.
      If there are missing values in the output due to math errors and 
      missing is set to None, output will fill those missing values 
      with the MA default value of 1e+20.

    Reference:
    * Wallace, J. M., and P. V. Hobbs (1977): Atmospheric Science:
      An Introductory Survey.  San Diego, CA:  Academic Press, ISBN
      0-12-732950-1, p. 52.

    Example without missing values:
    >>> from rho import rho
    >>> import Numeric as N
    >>> p  = N.array([1000.0, 800., 500.])
    >>> Tv = N.array([273.15, 280., 290.])
    >>> a = rho(p, Tv)
    >>> ['%.8g' % a[i] for i in range(len(a))]
    ['1.2753848', '0.99534675', '0.60064028']

    Example with missing values:
    >>> p  = N.array([1000.0, 1e+28, 500.])
    >>> Tv = N.array([273.15,  280., 290.])
    >>> a = rho(p, Tv, missing=1e+28)
    >>> ['%.8g' % a[i] for i in range(len(a))]
    ['1.2753848', '1e+28', '0.60064028']
    """
    import MA
    import Numeric as N
    from atmconst import AtmConst
    from is_numeric_float import is_numeric_float


    #- Check input is of the correct type:

    if is_numeric_float(p_in, Tv_in) != 1:
        raise TypeError, "rho:  Arg not Numeric floating"


    #- Change input to MA:  just set to input value unless there
    #  are missing values, in which case add mask:

    if missing == None:
        p  = MA.masked_array(p_in)
        Tv = MA.masked_array(Tv_in)
    else:
        p  = MA.masked_values(p_in, missing, copy=0)
        Tv = MA.masked_values(Tv_in, missing, copy=0)


    #- Load atmospheric constants:

    const = AtmConst()


    #- Calculate/return air density (eqn. 2.16, Wallace and Hobbs):

    return MA.filled( p * (100.0 / const.R_d) / Tv, missing )




#-------------------------- Main:  Test Module -------------------------

#- Define additional examples for doctest to use:

__test__ = { 'Additional Example 1':
    """
    >>> from rho import rho
    >>> import Numeric as N
    >>> p  = N.array([1000.0, 800., 500.])
    >>> Tv = N.array([273.15, 280., 1e+28])
    >>> a = rho(p, Tv, missing=1e+28)
    >>> ['%.8g' % a[i] for i in range(len(a))]
    ['1.2753848', '0.99534675', '1e+28']
    >>> p = N.array([850.0, 500.0, 250.0])
    >>> Tv = N.array([286.0715, 265.6786, 224.7649])
    >>> a = rho(p, Tv)
    >>> ['%.7g' % a[i] for i in range(len(a))]
    ['1.035111', '0.6556256', '0.3874842']
    """ }


#- Execute doctest if module is run from command line:

if __name__ == "__main__":
    """Test the module.

    Tests the examples in all the module documentation strings, plus 
    __test__.  The parent directory is added to sys.path for this 
    module testing case.
    """
    import doctest, sys, os
    sys.path.append(os.pardir)
    doctest.testmod(sys.modules[__name__])




# ===== end file =====