The fatiando package has been deprecated. Please check out the new tools in the Fatiando a Terra website: www.fatiando.org

# Source code for fatiando.gravmag.normal_gravity

"""
Gravity of ellipsoid models and common reductions (Bouguer, free-air)

**Reference ellipsoids**

This module uses instances of
:class:~fatiando.gravmag.normal_gravity.ReferenceEllipsoid to store the
physical constants of ellipsoids.
To create a new ellipsoid, just instantiate ReferenceEllipsoid and give it
the semimajor axis a, the flattening f, the geocentric gravitational
constant GM, and the angular velocity omega.
All  other quantities, like the gravity at the poles and equator,
eccentricities, etc, are computed by the class from these 4 parameters.

Available ellipsoids:

* WGS84 (values taken from Hofmann-Wellenhof and Moritz, 2005)::

>>> from fatiando.gravmag.normal_gravity import WGS84
>>> print(WGS84.name)
World Geodetic System 1984
>>> print('{:.0f}'.format(WGS84.a))
6378137
>>> print('{:.17f}'.format(WGS84.f))
0.00335281066474748
>>> print('{:.10g}'.format(WGS84.GM))
3.986004418e+14
>>> print('{:.7g}'.format(WGS84.omega))
7.292115e-05
>>> print('{:.4f}'.format(WGS84.b))
6356752.3142
>>> print('{:.8f}'.format(WGS84.E)) # Linear eccentricity
521854.00842339
>>> print('{:.15f}'.format(WGS84.e_prime)) # second eccentricity
0.082094437949696
>>> print('{:.10f}'.format(WGS84.gamma_a)) # gravity at the equator
9.7803253359
>>> print('{:.11f}'.format(WGS84.gamma_b)) # gravity at the pole
9.83218493786
>>> print('{:.14f}'.format(WGS84.m))
0.00344978650684

**Normal gravity**

* :func:~fatiando.gravmag.normal_gravity.gamma_somigliana: compute the normal
gravity using the Somigliana formula (Hofmann-Wellenhof and Moritz, 2005).
Calculated **on** the ellipsoid.
* :func:~fatiando.gravmag.normal_gravity.gamma_somigliana_free_air: compute
normal gravity at a height using the Somigliana formula plus the free-air
correction :math:-0.3086H\ mGal/m.
* :func:~fatiando.gravmag.normal_gravity.gamma_closed_form: compute normal
gravity using the closed form expression from Li and Gotze (2001). Can
compute anywhere (on, above, under the ellipsoid).

**Bouguer**

* :func:~fatiando.gravmag.normal_gravity.bouguer_plate: compute the
gravitational attraction of an infinite plate (Bouguer plate). Calculated
**on top** of the plate.

You can use :mod:fatiando.gravmag.prism and :mod:fatiando.gravmag.tesseroid
to calculate the terrain effect for a better correction.

**References**

Hofmann-Wellenhof, B. and H. Moritz, 2005, Physical Geodesy,
Springer-Verlag Wien, ISBN-13: 978-3-211-23584-3

Li, X. and H. J. Gotze, 2001, Tutorial: Ellipsoid, geoid, gravity, geodesy,
and geophysics, Geophysics, 66(6), p. 1660-1668, doi: 10.1190/1.1487109

----
"""
from __future__ import division
import math
import numpy

from .. import utils
from ..constants import G

[docs]class ReferenceEllipsoid(object):
"""
A generic reference ellipsoid.

It stores the physical constants defining the ellipsoid and has properties
for computing other (derived) quantities.

All quantities are expected and returned in SI units.

Parameters:

* a : float
The semimajor axis (the largest one, at the equator). In meters.
* f : float
* GM : float
The geocentric gravitational constant of the earth, including the
atmosphere. In :math:m^3 s^{-2}.
* omega : float
The angular velocity of the earth. In :math:rad s^{-1}.

"""

def __init__(self, name, a, f, GM, omega):
self.name = name
self._a = a
self._f = f
self._GM = GM
self._omega = omega

@property
def a(self):
"Semimajor axis"
return self._a

@property
def f(self):
"Flattening"
return self._f

@property
def GM(self):
"Geocentric gravitational constant (including the atmosphere)"
return self._GM

@property
def omega(self):
"Angular velocity"
return self._omega

@property
def b(self):
"Semiminor axis"
return self.a*(1 - self.f)

@property
def E(self):
"Linear eccentricity"
return math.sqrt(self.a**2 - self.b**2)

@property
def e_prime(self):
"Second eccentricity"
return self.E/self.b

@property
def m(self):
r":math:\omega^2 a^2 b / (GM)"
return (self.omega**2)*(self.a**2)*self.b/self.GM

@property
def gamma_a(self):
"Normal gravity at the equator"
bE = self.b/self.E
atanEb = math.atan2(self.E, self.b)
q0 = 0.5*((1 + 3*bE**2)*atanEb - 3*bE)
q0prime = 3*(1 + bE**2)*(1 - bE*atanEb) - 1
m = self.m
return self.GM*(1 - m - m*self.e_prime*q0prime/(6*q0))/(self.a*self.b)

@property
def gamma_b(self):
"Normal gravity at the poles"
bE = self.b/self.E
atanEb = math.atan2(self.E, self.b)
q0 = 0.5*((1 + 3*bE**2)*atanEb - 3*bE)
q0prime = 3*(1 + bE**2)*(1 - bE*atanEb) - 1
return self.GM*(1 + self.m*self.e_prime*q0prime/(3*q0))/self.a**2

WGS84 = ReferenceEllipsoid(a=6378137, f=1/298.257223563, GM=3986004.418e8,
omega=7292115e-11,
name="World Geodetic System 1984")

[docs]def gamma_somigliana(latitude, ellipsoid=WGS84):
'''
Calculate the normal gravity by using Somigliana's formula.

This formula computes normal gravity **on** the ellipsoid (height = 0).

Parameters:

* latitude : float or numpy array
The latitude where the normal gravity will be computed (in degrees)
* ellipsoid : :class:~fatiando.gravmag.normal_gravity.ReferenceEllipsoid
The reference ellipsoid used.

Returns:

* gamma : float or numpy array
The computed normal gravity (in mGal).

'''
sin2 = numpy.sin(lat)**2
cos2 = numpy.cos(lat)**2
top = ((ellipsoid.a*ellipsoid.gamma_a)*cos2
+ (ellipsoid.b*ellipsoid.gamma_b)*sin2)
bottom = numpy.sqrt(ellipsoid.a**2*cos2 + ellipsoid.b**2*sin2)
gamma = top/bottom
return utils.si2mgal(gamma)

[docs]def gamma_somigliana_free_air(latitude, height, ellipsoid=WGS84):
r'''
Calculate the normal gravity at a height using Somigliana's formula and the
free-air correction.

Parameters:

* latitude : float or numpy array
The latitude where the normal gravity will be computed (in degrees)
* height : float or numpy array
The height of computation (in meters). Should be ellipsoidal
(geometric) heights for geophysical purposes.
* ellipsoid : :class:~fatiando.gravmag.normal_gravity.ReferenceEllipsoid
The reference ellipsoid used.

Returns:

* gamma : float or numpy array
The computed normal gravity (in mGal).

'''
fa = -0.3086*height
gamma = gamma_somigliana(latitude, ellipsoid=ellipsoid) + fa
return gamma

[docs]def gamma_closed_form(latitude, height, ellipsoid=WGS84):
"""
Calculate normal gravity at a height using the closed form expression of
Li and Gotze (2001).

Parameters:

* latitude : float or numpy array
The latitude where the normal gravity will be computed (in degrees)
* height : float or numpy array
The height of computation (in meters). Should be ellipsoidal
(geometric) heights for geophysical purposes.
* ellipsoid : :class:~fatiando.gravmag.normal_gravity.ReferenceEllipsoid
The reference ellipsoid used.

Returns:

* gamma : float or numpy array
The computed normal gravity (in mGal).

"""
E2 = ellipsoid.E**2
bE = ellipsoid.b/ellipsoid.E
atanEb = math.atan2(ellipsoid.E, ellipsoid.b)
coslat = numpy.cos(lat)
sinlat = numpy.sin(lat)
tanlat = sinlat/coslat
beta = numpy.arctan2(ellipsoid.b*tanlat, ellipsoid.a)
sinbeta = numpy.sin(beta)
cosbeta = numpy.cos(beta)
zl2 = (ellipsoid.b*sinbeta + height*sinlat)**2
rl2 = (ellipsoid.a*cosbeta + height*coslat)**2
D = (rl2 - zl2)/E2
R = (rl2 + zl2)/E2
cosbetal = numpy.sqrt(0.5*(1 + R) - numpy.sqrt(0.25*(1 + R**2) - 0.5*D))
cosbetal2 = cosbetal**2
sinbetal2 = 1 - cosbetal2
bl = numpy.sqrt(rl2 + zl2 - E2*cosbetal2)
bl2 = bl**2
blE = bl/ellipsoid.E
atanEbl = numpy.arctan2(ellipsoid.E, bl)
q0 = 0.5*((1 + 3*bE**2)*atanEb - 3*bE)
q0l = 3*(1 + blE**2)*(1 - blE*atanEbl) - 1
W = numpy.sqrt((bl2 + E2*sinbetal2)/(bl2 + E2))
omega2 = ellipsoid.omega**2
a2 = ellipsoid.a**2
part1 = ellipsoid.GM/(bl2 + E2)
part2 = -cosbetal2*bl*ellipsoid.omega**2
part3 = 0.5*sinbetal2 - 1/6
part4 = a2*ellipsoid.E*q0l*omega2/((bl2 + E2)*q0)
gamma = utils.si2mgal((part1 + part2 + part3*part4)/W)
return gamma

[docs]def bouguer_plate(topography, density_rock=2670, density_water=1040):
r"""
Calculate the gravitational effect of an infinite Bouguer plate.

.. note:: The effect is calculated **on top** of the plate.

Uses the famous :math:g_{BG} = 2 \pi G \rho t formula, where t is the
height of the topography. On water (i.e., t < 0), uses
:math:g_{BG} = 2 \pi G (\rho_{water} - \rho_{rock})\times (-t).

Parameters:

* topography : float or numpy array
The height of topography (in meters).
* density_rock : float
The density of crustal rocks
* density_water : float
The density of ocean water

Returns:

* g_bouguer : float or array
The computed gravitational effect of the Bouguer plate

"""
t = numpy.atleast_1d(topography)
g_bg = numpy.empty_like(t)
g_bg[t >= 0] = 2*numpy.pi*G*density_rock*t[t >= 0]
g_bg[t < 0] = 2*numpy.pi*G*(density_water - density_rock)*(-t[t < 0])
g_bg = utils.si2mgal(g_bg)
if g_bg.size == 1:
return g_bg[0]
else:
return g_bg