add units and constants

docs/index.rst

@@ -8,6 +8,7 @@ Welcome to gMeterPy's documentation!
    :hidden:
 
    install.rst
+   user_guide.rst
    api/index.rst
    contributing.rst
 

docs/user_guide.rst

@@ -0,0 +1,71 @@
+User guide
+==========
+
+Units (`gmeterpy.units`)
+------------------------
+
+The unit aware calculations is one of the main goals of the `gMeterPy`. 
+We are using `astropy.units` sub-package for this. To use units you need 
+to import `gmeterpy.units` module:
+
+>>> import gmeterpy.units as u
+
+And now you can work with gravity acceleration SI units:
+
+>>> gravity = 9.8 * u.m / u.s**2
+>>> gravity
+<Quantity 9.8 m / s2>
+
+and convert them to more convinient CGS units:
+
+>>> gravity.to(u.Gal)
+<Quantity 980. Gal>
+
+and even use prefixes (u -- micro):
+
+>>> gravity.to(u.uGal)
+<Quantity 9.8e+08 uGal>
+
+In addition to Astropy built-in units we introduce Eotvos unit for
+the gravity gradient:
+
+>>> gradient = 0.3086 * u.mGal / u.m
+>>> gradient.to(u.E)
+<Quantity 3086. E>
+
+
+Constants (`gmeterpy.constants`)
+--------------------------------
+
+The `gMeterPy` uses `astropy.constants` sub-package for handling constants.
+We expand it for some frequently used constants in gravimetry.
+
++------------+-------+------------------------------------------------+
+| Constant   |  Description                                           |
++============+========================================================+
+| `omega`    | Mean angular rotation rate of the Earth                |
++------------+-------+------------------------------------------------+
+| `atm_sens` | Default gravity/pressure sensitivity (0.3 uGal / mbar) |
++------------+-------+------------------------------------------------+
+
+International Standard Atmosphere
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The Standard Atmosphere (ISO 2533:1975) is used to calculate 
+the normal pressure in the atmospheric pressure correction.
+
++----------+-------+----------------------------------------------+
+| Constant |  Description                                         |
++==========+======================================================+
+|   `p0`   | Standard atmospheric pressure at sea level           |
++----------+-------+----------------------------------------------+
+|   `L`    | Temperature lapse rate                               |
++----------+-------+----------------------------------------------+
+|   `Tn`   | Standard thermodynamic air temperature at sea level  |
++----------+-------+----------------------------------------------+
+|   `gn`   | Standard acceleration of free-fall                   |
++----------+-------+----------------------------------------------+
+|   `M`    | Air molar mass at sea level                          |
++----------+-------+----------------------------------------------+
+|   `R`    |  Universal gas constant                              |
++----------+-------+----------------------------------------------+

gmeterpy/constants.py

@@ -0,0 +1,50 @@
+# -*- coding: utf-8 -*-
+"""Constants for gMeterPy.
+
+"""
+
+from astropy.constants import Constant
+
+
+omega = Constant(
+    abbrev='omega',
+    name='Mean angular rotation rate of the Earth',
+    value=7292115.0e-11,
+    unit='rad / s',
+    uncertainty=0.0,
+    reference='')
+
+atm_sens = Constant(
+    abbrev='atm_sens',
+    name='Default sensitivity of changes in gravity '
+    'with variations in atmospheric pressure.',
+    value=0.3,
+    unit='uGal / mbar',
+    uncertainty=0.0,
+    reference='IAG 1983 Resolution no.9')
+
+
+# ISO 2533:1975 constants
+class ISA1975(Constant):
+    default_reference = 'Standard Atmosphere ISO 2533:1975'
+    _registry = {}
+    _has_incompatible_units = set()
+
+
+p0 = ISA1975('p0', 'Standard atmospheric pressure at sea level.',
+             101325.0, 'Pa', 0.0, system='si')
+
+L = ISA1975('L', 'Temperature lapse rate.',
+            0.0065, 'K / m', 0.0, system='si')
+
+Tn = ISA1975('Tn', 'Standard thermodynamic air temperature at sea level.',
+             288.15, 'K', 0.0, system='si')
+
+gn = ISA1975('gn', 'Standard acceleration of free-fall.',
+             9.80665, 'm / s**2', 0.0, system='si')
+
+M = ISA1975('M', 'Air molar mass at sea level.',
+            0.028964420, 'kg / mol', 0.0, system='si')
+
+R = ISA1975('R', 'Universal gas constant.',
+            8.31432, 'J / (mol * K)', 0.0, system='si')

gmeterpy/corrections/atmosphere.py

@@ -2,20 +2,26 @@
 """Atmospheric pressure correction.
 
 This module contains the atmospheric correction to the gravity observations.
+
 """
 
+import gmeterpy.units as u
+import gmeterpy.constants as c
+
 
-def normal_pressure(height):
-    r"""Normal atmospheric pressure, in Pa.
+@u.quantity_input
+def normal_pressure(height: u.m) -> u.Pa:
+    r"""Normal atmospheric pressure.
 
     Parameters
     ----------
-    height : float or array_like
-        Height above sea level, in metres.
+    height : ~astropy.units.Quantity
+        Height above sea level.
 
     Returns
     -------
-    float or array_like: Normal pressure, in Pa.
+    pn : ~astropy.units.Quantity
+        Normal pressure.
 
     Notes
     -----
@@ -24,7 +30,7 @@ def normal_pressure(height):
 
     .. math::
 
-       p_n = 1.01325\times 10^5 \left(1 -
+       p_n = 101325 \left(1 -
        0.0065 \dfrac{H}{288.15}\right)^{5.2559}\quad [\textrm{Pa}]
 
     where :math:`H` -- physical height of the station in metres.
@@ -35,29 +41,32 @@ def normal_pressure(height):
        1975 Standard Atmosphere ISO 2533 1975
 
     """
+    # pn = 101325 * (1 - 0.0065 * height / 288.15) ** 5.2559
+    pn = c.p0 * (1 - c.L * height / c.Tn)**(c.gn * c.M / c.R / c.L).round(4)
 
-    p_n = 101325 * (1 - 0.0065 * height / 288.15) ** 5.2559
-
-    return p_n
+    return pn
 
 
-def atmospheric_pressure_correction(height, pressure, barometric_factor=0.30):
-    r"""Atmospheric pressure correction, in m/s**2.
+@u.quantity_input
+def atmospheric_pressure_correction(height: u.m, pressure: u.Pa,
+                                    barometric_factor:
+                                    u.uGal / u.mbar = c.atm_sens) -> u.uGal:
+    r"""Atmospheric pressure correction.
 
     Parameters
     ----------
-    height : float or array_like
-        Height above sea level, in metres.
-    pressure : float or array_like
-        Observed atmospheric pressure, in Pa.
-    barometric_factor : float
-        Barometric factor, in 1e-10 m/s**2 / Pa.
-        Default is 0.30e-10 m/s**2 / Pa (or 0.3 muGal / mBar)
-        as recommended by IAG.
+    height : ~astropy.units.Quantity
+        Height above sea level.
+    pressure : ~astropy.units.Quantity
+        Observed atmospheric pressure.
+    barometric_factor : ~astropy.units.Quantity
+        Barometric factor.
+        Default is 0.3 uGal / mBar as recommended by IAG.
 
     Returns
     -------
-    float or array_like:  Atmospheric pressure correction, in m/s**2.
+    delta_g_atm : ~astropy.units.Quantity
+        Atmospheric pressure correction.
 
     Notes
     -----
@@ -68,14 +77,16 @@ def atmospheric_pressure_correction(height, pressure, barometric_factor=0.30):
 
     .. math::
 
-       \delta g = 0.30\times 10^{-10} (p_a - p_n)\quad [\textrm{ms}^{-2}]
+       \Delta g = 0.30\times 10^{-10} (p_a - p_n)\quad [\textrm{ms}^{-2}]
 
     where :math:`p_a` -- actual observed air pressure,
     :math:`p_n` -- normal pressure calculated by
-    :meth:`~gmeterpy.corrections.atmosphere.normal_pressure` function,
-    0.30 -- default barometric factor.
+    `gmeterpy.corrections.atmosphere.normal_pressure` function,
+    0.3 uGal / mBar -- default barometric factor.
 
-    """
+    The actual value of :math:`p_n` needs to be documented.
 
-    p_n = normal_pressure(height)
-    return barometric_factor * 1e-10 * (pressure - p_n)
+    """
+    delta_p = pressure - normal_pressure(height)
+    delta_g_atm = barometric_factor * delta_p
+    return delta_g_atm