Merge pull request #337 from nawendt/master

Added virtual temperature, virtual potential temperature, and density calculations.

docs/api/thermo.rst

@@ -21,3 +21,6 @@ Thermodynamic Calculations
     mixing_ratio
     saturation_mixing_ratio
     equivalent_potential_temperature
+    virtual_temperature
+    virtual_potential_temperature
+    density

metpy/calc/tests/test_thermo.py

@@ -6,10 +6,11 @@
 import numpy as np
 import pytest
 
-from metpy.calc import (dewpoint, dewpoint_rh, dry_lapse, equivalent_potential_temperature,
-                        lcl, lfc, mixing_ratio, moist_lapse, parcel_profile,
-                        potential_temperature, saturation_mixing_ratio,
-                        saturation_vapor_pressure, vapor_pressure)
+from metpy.calc import (density, dewpoint, dewpoint_rh, dry_lapse,
+                        equivalent_potential_temperature, lcl, lfc, mixing_ratio, moist_lapse,
+                        parcel_profile, potential_temperature, saturation_mixing_ratio,
+                        saturation_vapor_pressure, vapor_pressure,
+                        virtual_potential_temperature, virtual_temperature)
 from metpy.testing import assert_almost_equal, assert_array_almost_equal
 from metpy.units import units
 
@@ -190,3 +191,29 @@ def test_equivalent_potential_temperature():
     t = 288. * units.kelvin
     ept = equivalent_potential_temperature(p, t)
     assert_almost_equal(ept, 315.9548 * units.kelvin, 3)
+
+
+def test_virtual_temperature():
+    """Test virtual temperature calculation."""
+    t = 288. * units.kelvin
+    qv = .0016  # kg/kg
+    tv = virtual_temperature(t, qv)
+    assert_almost_equal(tv, 288.2796 * units.kelvin, 3)
+
+
+def test_virtual_potential_temperature():
+    """Test virtual potential temperature calculation."""
+    p = 999. * units.mbar
+    t = 288. * units.kelvin
+    qv = .0016  # kg/kg
+    theta_v = virtual_potential_temperature(p, t, qv)
+    assert_almost_equal(theta_v, 288.3620 * units.kelvin, 3)
+
+
+def test_density():
+    """Test density calculation."""
+    p = 999. * units.mbar
+    t = 288. * units.kelvin
+    qv = .0016  # kg/kg
+    rho = density(p, t, qv).to(units.kilogram / units.meter ** 3)
+    assert_almost_equal(rho, 1.2072 * (units.kilogram / units.meter ** 3), 3)

metpy/calc/thermo.py

@@ -513,3 +513,116 @@ def equivalent_potential_temperature(pressure, temperature):
     pottemp = potential_temperature(pressure, temperature)
     smixr = saturation_mixing_ratio(pressure, temperature)
     return pottemp * np.exp(Lv * smixr / (Cp_d * temperature))
+
+
+@exporter.export
+def virtual_temperature(temperature, mixing, molecular_weight_ratio=epsilon):
+    r"""Calculate virtual temperature.
+
+    This calculation must be given an air parcel's temperature and mixing ratio.
+    The implementation uses the formula outlined in [8]_.
+
+    Parameters
+    ----------
+    temperature: `pint.Quantity`
+        The temperature
+    mixing : `pint.Quantity`
+        dimensionless mass mixing ratio
+    molecular_weight_ratio : `pint.Quantity` or float, optional
+        The ratio of the molecular weight of the constituent gas to that assumed
+        for air. Defaults to the ratio for water vapor to dry air
+        (:math:`\epsilon\approx0.622`).
+
+    Returns
+    -------
+    `pint.Quantity`
+        The corresponding virtual temperature of the parcel
+
+    Notes
+    -----
+    .. math:: T_v = T \frac{\text{w} + \epsilon}{\epsilon\,(1 + \text{w})}
+
+    References
+    ----------
+    .. [8] Hobbs, Peter V. and Wallace, John M., 2006: Atmospheric Science, an Introductory
+           Survey. 2nd ed. 80.
+    """
+    return temperature * ((mixing + molecular_weight_ratio) /
+                          (molecular_weight_ratio * (1 + mixing)))
+
+
+@exporter.export
+def virtual_potential_temperature(pressure, temperature, mixing,
+                                  molecular_weight_ratio=epsilon):
+    r"""Calculate virtual potential temperature.
+
+    This calculation must be given an air parcel's pressure, temperature, and mixing ratio.
+    The implementation uses the formula outlined in [9]_.
+
+    Parameters
+    ----------
+    pressure: `pint.Quantity`
+        Total atmospheric pressure
+    temperature: `pint.Quantity`
+        The temperature
+    mixing : `pint.Quantity`
+        dimensionless mass mixing ratio
+    molecular_weight_ratio : `pint.Quantity` or float, optional
+        The ratio of the molecular weight of the constituent gas to that assumed
+        for air. Defaults to the ratio for water vapor to dry air
+        (:math:`\epsilon\approx0.622`).
+
+    Returns
+    -------
+    `pint.Quantity`
+        The corresponding virtual potential temperature of the parcel
+
+    Notes
+    -----
+    .. math:: \Theta_v = \Theta \frac{\text{w} + \epsilon}{\epsilon\,(1 + \text{w})}
+
+    References
+    ----------
+    .. [9] Markowski, Paul and Richardson, Yvette, 2010: Mesoscale Meteorology in the
+           Midlatitudes. 13.
+    """
+    pottemp = potential_temperature(pressure, temperature)
+    return virtual_temperature(pottemp, mixing, molecular_weight_ratio)
+
+
+@exporter.export
+def density(pressure, temperature, mixing, molecular_weight_ratio=epsilon):
+    r"""Calculate density.
+
+    This calculation must be given an air parcel's pressure, temperature, and mixing ratio.
+    The implementation uses the formula outlined in [10]_.
+
+    Parameters
+    ----------
+    temperature: `pint.Quantity`
+        The temperature
+    pressure: `pint.Quantity`
+        Total atmospheric pressure
+    mixing : `pint.Quantity`
+        dimensionless mass mixing ratio
+    molecular_weight_ratio : `pint.Quantity` or float, optional
+        The ratio of the molecular weight of the constituent gas to that assumed
+        for air. Defaults to the ratio for water vapor to dry air
+        (:math:`\epsilon\approx0.622`).
+
+    Returns
+    -------
+    `pint.Quantity`
+        The corresponding density of the parcel
+
+    Notes
+    -----
+    .. math:: \rho = \frac{p}{R_dT_v}
+
+    References
+    ----------
+    .. [10] Hobbs, Peter V. and Wallace, John M., 2006: Atmospheric Science, an Introductory
+           Survey. 2nd ed. 67.
+    """
+    virttemp = virtual_temperature(temperature, mixing, molecular_weight_ratio)
+    return (pressure / (Rd * virttemp)).to(units.kilogram / units.meter ** 3)