# World Geodetic System (1984)¶

The World Geodetic System 1984 (WGS 84) [WGS84] describes the best geodetic reference system for the Earth available for the practical applications of mapping, charting, geopositioning, and navigation, using data, techniques and technology available through 2013 by the United States of America’s National Geospatial-Intelligence Agency (NGA.)

The WGS 84 Coordinate System is a Conventional Terrestrial Reference System (CTRS), that follows the criteria outlined in the International Earth Rotation and Reference Systems Service (IERS):

• It is geocentric, meaning the center of mass is defined for the whole Earth including oceans and atmosphere.
• Its scale is that of the local Earth frame, in the meaning of a relativistic theory of gravitation.
• Its orientation was initially given by the Bureau International de l’Heure (BIH) at the epoch 1984.0
• It is defined as a right-handed, orthogonal, and Earth-fixed coordinate system, where the Z-axis serves as the rotational axis of the ellipsoid revolution.

WGS 84 (G1762) is the sixth and latest update (16 Oct 2013) to the realization of the WGS 84 Reference Frame, and is updated to incorporate international conventions and alignements to the International Terrestrial Reference Frame 2008 (ITRF2008.)

A simple and practical approach is a mathematically manageable reference surface (an ellipsoid) approximating the broad features of the figure and of the gravity field of the Earth, which is used to describe the Earth surface and the forces acting on and above it.

A second, and more precise, solution considers an equipotential surface of the gravity field of the Earth, called the geoid, which is used to define the Earth Gravitational Model (EGM.) The implementation of this second model falls out of the scope of this module.

The WGS 84 ellipsoid is defined by the semi-major axis ($$a$$), the reciprocal flattening ($$1/f$$) of an oblate ellipsoid of revolution, the Geocentric Gravitational Constant ($$GM$$) and the angular velocity ($$\omega$$)

Symbol Definition Value Unit
$$a$$ Semi-major Axis 6378137.0 m
$$1/f$$ Flattening Factor of the Earth 298.257223563
$$GM$$ Geocentric Gravitational Constant 3.986004418 x 10^14 m^3/s^2
$$\omega$$ Earth’s Nominal Mean Angular Velocity 7.292115 x 10^-5 rad/s

The first two parameters ($$a$$, $$1/f$$) define the geometry of the rotational ellipsoid, while the other two parameters ($$GM$$, $$\omega$$) permit the unique determination of its associated normal gravity field.

Having these 4 elements defined, it is possible to estimate most of the WGS84 parameters directly. The class WGS already sets these values by default (and estimates the semi-minor axis b), but with slightly different notation:

>>> wgs = ahrs.utils.WGS()
>>> wgs.a       # Semi-major Axis
6378137.0
>>> 1/wgs.f     # Flattening Factor of the Earth
298.257223563
>>> wgs.gm      # Geocentric Gravitational Constant
398600441800000.0
>>> wgs.w       # Earth's Nominal Mean Angular Velocity
7.292115e-05
>>> wgs.b       # Semi-minor Axis
6356752.314245179


Furthermore, there are three unitless values used in the computation of the Moments of Inertia:

Symbol Definition Value
C2,0 dyn[2008] Dynamic Second Degree Zonal Harmonics -4.84165143790815 x 10^-4
C2,2 dyn[2008] Dynamic Second Degree Sectorial Harmonics 2.43938357328313 x 10^-6
H Dynamic Ellipticity 3.2737949 x 10^-3

The dynamic harmonics are recovered from the data obtained empirically in EGM2008, and are NOT derived from the ellipsoid parameters. The subscript dyn[2008] denote their origin.

The dynamic ellipticity is a factor in the theoretical value of the rate of precession of the equinoxes, which is known from observation too.

The WGS 84 Ellipsoid is identified as a geocentric, equipotential ellipsoid of revolution, i.e., an ellipsoid with a surface on which the value of the gravity potential is the same everywhere.

## Earth’s Gravity Field¶

In a rectangular system a point $$\mathbf{p}$$ is located with the coordinates $$\begin{pmatrix}x & y & z\end{pmatrix}$$, but in ellipsoidal coordinates, this point is found with $$\begin{pmatrix}u & \phi & \lambda\end{pmatrix}$$, where $$u$$ is the semi-minor axis of the ellipsoid, $$\phi$$ is the angle between the plumb line and the equatorial plane called the geographical latitude, and $$\lambda$$ is the angle between $$\mathbf{p}$$ and the meridian plane of Greenwich called the geographical longitude.

We assume Earth is an ellipsoid of revolution which is an equipotential surface of a normal gravity field. Although the Earth is not a perfect ellipsoid, its gravity field is easier to handle assuming it is one. The deviations of the field are so small that they are considered linear for this model.

The definition of the potential of the normal gravity field $$U$$ [Heiskanen] is:

$U = V + \Phi$

where $$V$$ is the potential of gravitational force defined as [1]:

$V = \frac{GM}{E}\arctan\Big(\frac{E}{u}\Big) + \frac{1}{2}\omega^2a^2\frac{q}{q_0}\Big(\sin^2\beta-\frac{1}{3}\Big)$

and $$\Phi$$ is the potential of centrifugal force:

$\Phi = \frac{1}{2} \omega^2(u^2+E^2)\cos^2\beta$

Here, the ellipsoid’s linear eccentricity is simply $$E = \sqrt{a^2-b^2}$$, and we set the helper constants:

$\begin{split}\begin{array}{ll} q &= \frac{1}{2} \Big[\Big(1+\frac{3u^2}{E^2}\Big)\arctan\frac{E}{u}-\frac{3u}{E}\Big] \\ q_0 &= \frac{1}{2} \Big[\Big(1+\frac{3b^2}{E^2}\Big)\arctan\frac{E}{b}-\frac{3b}{E}\Big] \\ q_0' &= 3\Big[\Big(1+\frac{1}{e'^2}\Big)\Big(1-\frac{1}{e'}\arctan e'\Big)\Big] - 1 \end{array}\end{split}$

The normal gravity vector $$\vec{\mathbf{g}}$$ is obtained from its potential $$U$$ :

$\begin{split}\vec{\mathbf{g}} = \nabla U = \begin{bmatrix}g_u \\ g_\beta \\ g_\lambda \end{bmatrix} = \begin{bmatrix} \frac{1}{w}\frac{\partial U}{\partial u} \\ \frac{1}{w\sqrt{u^2+E^2}}\frac{\partial U}{\partial \beta} \\ 0 \end{bmatrix}\end{split}$

We see that $$g_\lambda=0$$, because the effects of the potential along the longitudinal direction are neglected.

## Normal Gravity on the Surface¶

The magnitude of the normal gravity vector, simply called normal gravity, is:

$g = \|\vec{\mathbf{g}}\| = \sqrt{g_u^2+g_\beta^2}$

At the surface of the ellipsoid $$u=b$$ making $$g_{\beta,0}=0$$. So, the total gravity on the surface of the ellipsoid $$g_0$$ is just:

$g_0 = \|\vec{\mathbf{g}_0}\| = \frac{GM}{a\sqrt{a^2\sin^2\beta+b^2\cos^2\beta}} \Big[\Big(1+\frac{me'q_0'}{3q_0}\Big)\sin^2\beta + \Big(1-m-\frac{me'q_0'}{6q_0}\Big)\cos^2\beta\Big]$

where:

$m = \frac{\omega^2a^2b}{GM}$

and $$e'=\frac{E}{b}$$ is the second eccentrictiy.

From here we can find two basic variables that will help us to linearly estimate the normal gravity at any point on the ellipsoid.

Keeping it on the surface, we can estimate the normal gravity at the equator:

$g_e = g_0(\beta=0°) = \frac{GM}{ab}\Big(1-m-\frac{me'q_0'}{6q_0}\Big)$

Similarly, we estimate the normal gravity at the poles:

$g_p = g_0(\beta=90°) = \frac{GM}{a^2}\Big(1+\frac{me'q_0'}{3q_0}\Big)$

With these two basic values, we can linearly approximate the normal at any latitude on the surface of the ellipsoid, but we need to do it in geographical coordinates.

Using the property $$\tan\beta=\frac{b}{a}\tan\phi$$, we define a closed form formula to find the normal gravity $$g$$ at any given latitude $$\phi$$ [Somigliana1929]:

$g(\phi) = \frac{ag_e \cos^2\phi + bg_p\sin^2\phi}{\sqrt{a^2\cos^2\phi + b^2\sin^2\phi}}$

For numerical computation, a more convenient form is:

$g(\phi) = g_e\frac{1+k\sin^2\phi}{\sqrt{1-e^2\sin^2\phi}}$

using the helper variable $$k = \frac{bg_p}{ag_e} - 1$$.

The estimation of the normal gravity is already simplified with the class WGS requiring the latitude only.

>>> wgs.normal_gravity(0.0)     # Normal gravity at Equator (latitude = 0.0 °)
9.78032533590406
>>> wgs.normal_gravity(50.0)    # Normal gravity at latitude = 50.0 °
9.810702135603085


## Normal Gravity above the Surface¶

At small heights above the surface, the normal gravity can be approximated with a truncated Taylor Series with a positive direction downward along the geodetic normal to the ellipsoid:

$g(\phi, h) = g(\phi) \Big(1 - \frac{2}{a}\big(1+f+m-2f\sin^2\phi\big)h + \frac{3}{a^2}h^2\Big)$

where $$h$$ is the height, in meters, above the ellipsoid’s surface.

>>> wgs.normal_gravity(50.0, 1000.0)    # Gravity at latitude = 50.0 °, 1000 m above surface
9.807617683884756


## Other Gravitational Methods¶

The well known International Gravity Formula [Lambert] as described by Helmut Moritz in [Tscherning] for the Geodetic Reference System 1980 is also implemented here:

$g(\phi) = 9.780327 (1 + 0.0053024 \sin^2\phi - 0.0000058 \sin^2(2\phi))$
>>> ahrs.utils.international_gravity(10.0)
9.781884110728155


As a bonus, the normal gravity estimation of the European Cooperation on Legal Metrology (WELMEC) is also implemented here:

$g(\phi, h) = 9.780318(1 + 0.0053024\sin^2(\phi) - 0.0000058\sin^2(2\phi)) - 0.000003085h$
>>> ahrs.utils.welmec_gravity(50.0, 1000.0)     # 50.0° N, 1000 m above sea level
9.807610187885896


Although this is thought and mainly used for European latitudes.

All methods above can be used for cartography and basic gravitational references. For more advanced and precise estimations it is recommended to use the EGM2008, which is defined, like the WMM, with spherical harmonics, but of degree 2190 and order 2159.

Because of its high complexity and demanding computation, the implementation of the EGM2008 is left out of this module in favour of the more convenient applications developed in Fortran by the NGA.

## Footnotes¶

 [1] The most precise EGM2008 defines the potential of gravitational force in terms of spherical harmonics. But that definition is out of the scope of this package.

References

 [WGS84] World Geodetic System 1984. Its Definition and Relationships with Local Geodetic Systems. National Geospatial-Intelligence Agency (NGA) Standarization Document. 2014. (ftp://ftp.nga.mil/pub2/gandg/website/wgs84/NGA.STND.0036_1.0.0_WGS84.pdf)
 [Heiskanen] Heiskanen, W. A. and Moritz, H. Physical Geodesy. W. H. Freeman and Company. 1967.
 [WELMEC2009] WELMEC DIRECTIVE 2009/23/EC: Common application non-automatic weighing instruments. (https://www.welmec.org/documents/guides/2/)
 [Lambert] (1, 2) Walter D. Lambert. The International Gravity Formula. U.S. Coast and Geodetic Survey. 1945. (http://earth.geology.yale.edu/~ajs/1945A/360.pdf)
 [Somigliana1929] Carlo Somigliana. Teoria generale del campo gravitazionale dell’ellissoide di rotazione. Memorie della Società Astronomia Italiana, Vol. 4, p.425. 1929 (http://articles.adsabs.harvard.edu/pdf/1929MmSAI…4..425S)
 [Tscherning] (1, 2) C. Tscherning. The Geodesist’s Handbook 1984. Association Internationale de Géodésie. 1984. (https://office.iag-aig.org/doc/5d7f91ee333a3.pdf)
class ahrs.utils.wgs84.WGS(a: float = 6378137.0, f: float = 0.0033528106647474805, GM: float = 398600441800000.0, w: float = 7.292115e-05)

World Geodetic System 1984

Parameters: a (float, default: 6378137.0) – Objects’s Semi-major axis (Equatorial Radius), in meters. Defaults to Earth’s semi-major axis. f (float, default: 0.0033528106647474805) – Objects’s flattening factor. Defaults to Earth’s flattening. GM (float, default: 3.986004418e14) – Objects’s Standard Gravitational Constant in m^3/s^2 w (float, default: 0.00007292115) – Objects’s rotation rate in rad/s a (float) – Ellipsoid’s semi-major axis (Equatorial Radius), in meters. f (float) – Ellipsoid’s flattening factor. gm (float) – Objects’s Standard Gravitational Constant in m^3/s^2. w (float) – Ellipsoid’s rotation rate in rad/s. b (float) – Ellipsoid’s semi-minor axis (Polar Radius), in meters. is_geodetic (bool) – Whether the object describes Earth.
arithmetic_mean_radius
$R_1 = a\Big(1-\frac{f}{3}\Big)$

Example

>>> wgs = ahrs.utils.WGS()
6371008.771415059

Type: Mean Radius $$R_1$$ of the Three Semi-Axes, computed as
aspect_ratio
$AR = \frac{b}{a}$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.aspect_ratio
0.9966471893352525

Type: Aspect Ratio $$AR$$, computed as
atmosphere_gravitational_constant
$GM_A = G\ M_A$

where $$M_A$$ is the total mean mass of the atmosphere (with water vapor) equal to $$5.148\times 10^{18} \mathrm{kg}$$.

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.atmosphere_gravitational_constant
343591934.4

Type: Gravitational Constant of the Atmosphere $$GM_A$$, computed as
authalic_sphere_radius
$R_2 = R_P \Big(1-\frac{2}{3}e'^2 + \frac{26}{45}e'^4 - \frac{100}{189}e'^6 + \frac{7034}{14175}e'^8 - \frac{220652}{467775}e'^{10} \Big)$

Example

>>> wgs = ahrs.utils.WGS()
6371007.1809182055

Type: Radius $$R_2$$ of a Sphere of Equal Area, computed as
curvature_polar_radius
$\begin{split}\begin{array}{ll} R_P &= \frac{a^2}{b} = \frac{a}{\sqrt{1-e^2}} \\ &= \frac{a}{1-f} \end{array}\end{split}$

Example

>>> wgs = ahrs.utils.WGS()
6399593.625758493

Type: Polar Radius of Curvature $$R_P$$, computed as
dynamic_inertial_moment_about_X

Dynamic Moment of Inertia ($$A$$), with respect to the X-Axis of Rotation, computed as:

$A_{dyn} = -\sqrt{5}Ma^2\Big[\Big(1-\frac{1}{H}\Big)C_{2,0dyn} - \frac{C_{2,2dyn}}{\sqrt{3}}\Big]$

where $$H=3.2737949 \times 10^{-3}$$ is the Dynamic Ellipticity.

Example

>>> wgs = ahrs.utils.WGS()
8.007921777277886e+37

dynamic_inertial_moment_about_Y

Dynamic Moment of Inertia ($$B$$), with respect to the X-Axis of Rotation, computed as:

$B_{dyn} = -\sqrt{5}Ma^2\Big[\Big(1-\frac{1}{H}\Big)C_{2,0dyn} + \frac{C_{2,2dyn}}{\sqrt{3}}\Big]$

where $$H=3.2737949 \times 10^{-3}$$ is the Dynamic Ellipticity.

Example

>>> wgs = ahrs.utils.WGS()
8.008074799852911e+37

dynamic_inertial_moment_about_Z

Dynamic Moment of Inertia ($$C$$), with respect to the Z-Axis of Rotation, computed as:

$C_{dyn} = -\sqrt{5}Ma^2\frac{C_{2,0dyn}}{H}$

where $$H=3.2737949 \times 10^{-3}$$ is the Dynamic Ellipticity.

Example

>>> wgs = ahrs.utils.WGS()
8.03430094201443e+37

dynamical_form_factor
$J_2 = \frac{e^2}{3} \Big(1-\frac{2me'}{15q_0}\Big)$

where:

$q_0 = \frac{1}{2}\Big[\Big(1+\frac{3}{e'^2}\Big)\arctan e' - \frac{3}{e'}\Big]$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.dynamical_form_factor
0.0010826298213129219

Type: WGS 84 Dynamical Form Factor $$J_2$$, computed as
equatorial_normal_gravity

Normal Gravity $$g_e$$ at the Equator, in $$\frac{\mathrm{m}}{\mathrm{s}^2}$$, computed as:

$g_e = \frac{GM}{ab}\Big(1-m-\frac{me'q_0'}{6q_0}\Big)$

where:

$\begin{split}\begin{array}{ll} q_0 &= \frac{1}{2}\Big[\Big(1+\frac{3}{e'^2}\Big)\arctan e' - \frac{3}{e'}\Big] \\ q_0' &= 3\Big[\Big(1+\frac{1}{e'^2}\Big)\Big(1-\frac{1}{e'}\arctan e'\Big)\Big] - 1 \end{array}\end{split}$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.equatorial_normal_gravity
9.78032533590406

equivolumetric_sphere_radius
$\begin{split}\begin{array}{ll} R_3 &= \sqrt[3]{a^2b} \\ &= a\sqrt[3]{1-f} \end{array}\end{split}$

Example

>>> wgs = ahrs.utils.WGS()
6371000.790009159

Type: Radius $$R_3$$ of a Sphere of Equal Volume, computed as
first_eccentricity_squared
$e^2 = 2f - f^2$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.first_eccentricity_squared
0.0066943799901413165

Type: First Eccentricity Squared $$e^2$$ of the ellipsoid, computed as
geometric_dynamic_ellipticity
$H_{geo} = \frac{C_{geo}-A_{geo}}{C_{geo}}$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.geometric_dynamic_ellipticity
0.003258100628533992

Type: Geometric Solution for Dynamic Ellipticity $$H$$, computed as
geometric_inertial_moment

Geometric Moment of Inertia ($$A$$), with respect to Any Axis in the Equatorial Plane, computed as:

$A_{geo} = C_{geo} + \sqrt{5}Ma^2 C_{2,0geo}$

where $$C_{2,0geo} = -4.84166774985\times 10^{-4}$$ is Earth’s Geographic Second Degree Zonal Harmonic.

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.geometric_inertial_moment
8.046726628049449e+37

geometric_inertial_moment_about_Z

Geometric Moment of Inertia ($$C$$), with respect to the Z-Axis of Rotation, computed as:

$C_{geo} = \frac{2}{3}Ma^2\Big(1-\frac{2}{5}\sqrt{\frac{5m}{2f}-1}\Big)$

Example

>>> wgs = ahrs.utils.WGS()
8.073029370114392e+37

gravitational_constant_without_atmosphere

Geocentric Gravitational Constant with Earth’s Atmosphere Excluded $$GM'$$, computed as:

$GM' = GM - GM_A$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.gravitational_constant_without_atmosphere
398600098208065.6

linear_eccentricity
$E = \sqrt{a^2-b^2}$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.linear_eccentricity
521854.00842338527

Type: Linear Eccentricity $$E$$, computed as
mass
$M = \frac{GM}{G}$

where $$G$$ is the universal constant of gravitation equal to $$6.67428\times 10^{-11} \frac{\mathrm{m}^3}{\mathrm{kg}\ \mathrm{s}^2}$$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.mass
5.972186390142457e+24

Type: The Mass $$M$$ of the Earth, in kg, computed as
mean_normal_gravity

Mean Value $$\bar{g}$$ of Normal Gravity, in $$\frac{\mathrm{m}}{\mathrm{s}^2}$$, computed as:

$\bar{g} = g_e\Big(1 + \frac{1}{6}e^2 + \frac{1}{3}k + \frac{59}{360}e^4 + \frac{5}{18}e^2k + \frac{2371}{15120}e^6 + \frac{259}{1080}e^4k + \frac{270229}{1814400}e^8 + \frac{9623}{45360}e^6k \Big)$

where:

$k = \frac{bg_p - ag_e}{ag_e} = \frac{bg_p}{ag_e}-1$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.mean_normal_gravity
9.797643222256516

meridian_curvature_radius(lat: float) → float

Radius of the curvature in the prime meridian, computed as:

$R_M = \frac{a(1-e^2)}{\sqrt[3]{1-e^2\sin^2\phi}}$
Parameters: lat (float) – Geographical latitude, in decimal degrees.
normal_gravity(lat: float, h: float = 0.0) → float

Normal Gravity on (or above) Ellipsoidal Surface

Estimate the normal gravity on or above the surface of an ellipsoidal body using Somigliana’s formula (on surface) and a series expansion (above surface).

Somigliana’s closed formula as desribed by H. Moritz in [Tscherning] is:

$g = \frac{ag_e \cos^2\phi + bg_p\sin^2\phi}{\sqrt{a^2cos^2\phi + b^2\sin^2\phi}}$

For numerical computation, a more convenient form is:

$g = g_e\frac{1+k\sin^2\phi}{\sqrt{1-e^2\sin^2\phi}}$

with the helper constant $$k$$:

$k = \frac{bg_p}{ag_e}-1$
Parameters: lat (float) – Geographical latitude, in decimal degrees. h (float, default: 0.0) – Mean sea level height, in meters. g – Normal gravity at given point in space, in m/s^2. float

Examples

>>> wgs = ahrs.utils.WGS()
>>> wgs.normal_gravity(50.0)
9.810702135603085
>>> wgs.normal_gravity(50.0, 100.0)
9.810393625316983

normal_gravity_constant
$m = \frac{\omega^2a^2b}{GM}$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.normal_gravity_constant
0.0034497865068408447

Type: Normal Gravity Formula Constant $$m$$, computed as
normal_gravity_potential

Normal Gravity Potential $$U_0$$ of the WGS 84 Ellipsoid, computed as:

$U_0 = \frac{GM}{E} \arctan e' + \frac{\omega^2a^2}{3}$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.normal_gravity_potential
62636851.71456948

polar_normal_gravity

Normal Gravity $$g_p$$ at the Pole, in $$\frac{\mathrm{m}}{\mathrm{s}^2}$$, computed as:

$g_p = \frac{GM}{a^2}\Big(1+\frac{me'q_0'}{3q_0}\Big)$

where:

$\begin{split}\begin{array}{ll} q_0 &= \frac{1}{2}\Big[\Big(1+\frac{3}{e'^2}\Big)\arctan e' - \frac{3}{e'}\Big] \\ q_0' &= 3\Big[\Big(1+\frac{1}{e'^2}\Big)\Big(1-\frac{1}{e'}\arctan e'\Big)\Big] - 1 \end{array}\end{split}$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.polar_normal_gravity
9.832184937863065

second_degree_zonal_harmonic
$C_{2,0} = -\frac{J_2}{\sqrt{5}}$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.second_degree_zonal_harmonic
-0.00048416677498482876

Type: WGS 84 Second Degree Zonal Harmonic $$C_{2,0}$$, computed as
second_eccentricity_squared
$e'^2 = \frac{a^2-b^2}{b^2}$

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.second_eccentricity_squared
0.006739496742276434

Type: Second Eccentricity Squared $$e'^2$$, computed as
vertical_curvature_radius(lat: float) → float

Radius of the curvature in the prime vertical, computed as:

$R_N = \frac{a}{\sqrt{1-e^2\sin^2\phi}}$
Parameters: lat (float) – Geographical latitude, in decimal degrees.
ahrs.utils.wgs84.international_gravity(lat: float, epoch: str = '1980') → float

International Gravity Formula

Estimate normal gravity with the International Gravity Formula [Lambert], adapted from Stokes’ formula, and adopted by the International Association of Geodesy at its Stockholm Assembly in 1930.

The expression for gravity on a spheroid, which combines gravitational attraction and centrifugal acceleration, can be written in the form of a series:

$g = g_e\big(1 + \beta\sin^2(\phi) - \beta_1\sin^2(2\phi) - \beta_2\sin^2(\phi)\sin^2(2\phi) - \beta_3\sin^4(\phi)\sin^2(2\phi) - \dots\big)$

where the values of the $$\beta$$’s are:

$\begin{split}\begin{array}{ll} \beta &= \frac{5}{2}m\Big(1-\frac{17}{35}f - \frac{1}{245}f^2 - \frac{13}{18865}f^3 - \dots\Big) - f \\ \beta_1 &= \frac{1}{8}f(f+2\beta) \\ \beta_2 &= \frac{1}{8}f^2(2f+3\beta) - \frac{1}{32}f^3(3f+4\beta) \\ & \vdots \\ & \mathrm{etc.} \end{array}\end{split}$

For the case of the International Ellipsoid, the third-order terms are negligible. So, in practice, the term $$\beta_2$$ and all following terms are dropped to yield the form:

$g = g_e \big(1 + \beta \sin^2\phi - \beta_1 \sin^2(2\phi)\big)$

In the original definition the values of $$\beta$$ and $$\beta_1$$ are rounded off to seven decimal places to simply get the working formula:

$g = 9.78049 \big(1 + 0.0052884 \sin^2\phi - 0.0000059 \sin^2(2\phi)\big)$

Originally, the definitions of the elementary properties ($$a$$, $$g_e$$, etc.) weren’t as accurate as now. At different moments in history, the values were updated to improve the accuracy of the formula. Those different moments are named epochs and are labeled according to the year they were updated:

epoch $$g_e$$ $$\beta$$ $$\beta_1$$
1930 9.78049 5.2884 x 10^-3 5.9 x 10^-6
1948 9.780373 5.2891 x 10^-3 5.9 x 10^-6
1967 9.780318 5.3024 x 10^-3 5.9 x 10^-6
1980 9.780327 5.3024 x 10^-3 5.8 x 10^-6

The latest epoch, 1980, is used here by default.

Parameters: lat (float) – Geographical Latitude, in decimal degrees. epoch (str, default: '1980') – Epoch of the Geodetic Reference System. Options are '1930', '1948', '1967' and '1980'. g – Normal gravity, in m/s^2, at given latitude. float

Examples

>>> ahrs.utils.international_gravity(10.0)
9.781884110728155
>>> ahrs.utils.international_gravity(10.0, epoch='1930')
9.7820428934191

ahrs.utils.wgs84.welmec_gravity(lat: float, h: float = 0.0) → float

Reference normal gravity of WELMEC’s gravity zone

Gravity zones are implemented by European States on their territories for weighing instruments that are sensitive to gravity variations [WELMEC2009].

Manufacturers may adjust their instruments using the reference gravity formula:

$g = 9.780318(1 + 0.0053024\sin^2(\phi) - 0.0000058\sin^2(2\phi)) - 0.000003085h$

where $$\phi$$ is the geographical latitude and $$h$$ is the height in meters.

Parameters: lat (float) – Geographical Latitude, in decimal degrees. h (float, default: 0.0) – Mean sea level height, in meters. g – Normal gravity at given point in space, in m/s^2. float

Examples

>>> ahrs.utils.welmec_gravity(52.3, 80.0)      # latitude = 52.3°, height = 80 m
9.818628439187075