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
C 2,0 dyn[2008]	Dynamic Second Degree Zonal Harmonics	-4.84165143790815 x 10^-4
C 2,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}\mathrm{g}_u \\ \mathrm{g}_\beta \\ \mathrm{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 \(\mathrm{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:

\[\mathrm{g} = \|\vec{\mathbf{g}}\| = \sqrt{\mathrm{g}_u^2+\mathrm{g}_\beta^2}\]

where \(u=b\) at the surface of the ellipsoid making \(\mathrm{g}_{\beta,0}=0\). So, the total gravity on its surface, \(\mathrm{g}_0\), is just:

\[\mathrm{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 eccentricity.

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:

\[\mathrm{g}_e = \mathrm{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:

\[\mathrm{g}_p = \mathrm{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 \(\mathrm{g}\) at any given latitude \(\phi\) [Somigliana1929]:

\[\mathrm{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:

\[\mathrm{g}(\phi) = \mathrm{g}_e\frac{1+k\sin^2\phi}{\sqrt{1-e^2\sin^2\phi}}\]

using the helper variable \(k = \frac{b\mathrm{g}_p}{a\mathrm{g}_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:

\[\mathrm{g}(\phi, h) = \mathrm{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 implemented in ahrs:

\[\mathrm{g}(\phi) = 9.780327 (1 + 0.0053024 \sin^2\phi - 0.0000058 \sin^2(2\phi))\]

>>> ahrs.utils.international_gravity(10.0)
9.781884110728155

Additionally, the normal gravity estimation of the European Cooperation on Legal Metrology (WELMEC) is also included here:

\[\mathrm{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.

[WELMEC2021]

WELMEC Directives 2014/31/EU and 2014/32/EU: Common Application. Non-Automatic Weighing Instruments (NAWI); Automatic Weighing Instruments (AWI); Multi-dimensional Measuring Instruments (MDMI.) (https://www.welmec.org/welmec/documents/guides/2/2021/WELMEC_Guide_2_v2021.pdf)

[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) – Ellipsoid’s Semi-major axis (Equatorial Radius), in meters. Defaults to Earth’s semi-major axis.

• f (float, default: 0.0033528106647474805) – Ellipsoid’s flattening factor. Defaults to Earth’s flattening.

• GM (float, default: 3.986004418e14) – Ellipsoid’s Standard Gravitational Constant in m^3/s^2

• w (float, default: 0.00007292115) – Ellipsoid’s rotation rate in rad/s

property arithmetic_mean_radius: float

Mean Radius \(R_1\) of the Three Semi-Axes, computed as:

\[R_1 = a\Big(1-\frac{f}{3}\Big)\]

Example

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

property aspect_ratio: float

Aspect Ratio \(AR\), computed as:

\[AR = \frac{b}{a}\]

Example

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

property atmosphere_gravitational_constant: float

Gravitational Constant of the Atmosphere \(GM_A\), computed as:

\[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

property authalic_sphere_radius: float

Radius \(R_2\) of a Sphere of Equal Area, computed as:

\[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()
>>> wgs.authalic_sphere_radius
6371007.1809182055

property curvature_polar_radius: float

Polar Radius of Curvature \(R_P\), computed as:

\[\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()
>>> wgs.curvature_polar_radius
6399593.625758493

property dynamic_inertial_moment_about_X: float

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()
>>> wgs.dynamic_inertial_moment_about_X
8.007921777277886e+37

property dynamic_inertial_moment_about_Y: float

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()
>>> wgs.dynamic_inertial_moment_about_Y
8.008074799852911e+37

property dynamic_inertial_moment_about_Z: float

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()
>>> wgs.dynamic_inertial_moment_about_Z
8.03430094201443e+37

property dynamical_form_factor: float

WGS 84 Dynamical Form Factor \(J_2\), computed as:

\[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

property equatorial_normal_gravity: float

Normal Gravity \(\mathrm{g}_e\) at the Equator, in \(\frac{\mathrm{m}}{\mathrm{s}^2}\), computed as:

\[\mathrm{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

property equivolumetric_sphere_radius: float

Radius \(R_3\) of a Sphere of Equal Volume, computed as:

\[\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()
>>> wgs.equivolumetric_sphere_radius
6371000.790009159

property first_eccentricity_squared: float

First Eccentricity Squared \(e^2\) of the ellipsoid, computed as:

\[e^2 = 2f - f^2\]

Example

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

property geometric_dynamic_ellipticity: float

Geometric Solution for Dynamic Ellipticity \(H\), computed as:

\[H_{geo} = \frac{C_{geo}-A_{geo}}{C_{geo}}\]

Example

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

property geometric_inertial_moment: float

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

property geometric_inertial_moment_about_Z: float

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()
>>> wgs.geometric_inertial_moment_about_Z
8.073029370114392e+37

property gravitational_constant_without_atmosphere: float

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

property is_geodetic: bool

Check whether the ellipsoid model is geodetic (describes planet Earth.)

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.is_geodetic
True
>>> wgs = ahrs.utils.WGS(a=6_500_000)
>>> wgs.is_geodetic
False

property linear_eccentricity: float

Linear Eccentricity \(E\), computed as:

\[E = \sqrt{a^2-b^2}\]

Example

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

property mass: float

The Mass \(M\) of the Earth, in kg, computed as:

\[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

property mean_normal_gravity: float

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

\[\bar{\mathrm{g}} = \mathrm{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{b\mathrm{g}_p - a\mathrm{g}_e}{a\mathrm{g}_e} = \frac{b\mathrm{g}_p}{a\mathrm{g}_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, estimated at a given latitude, \(\phi\), 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:

\[\mathrm{g} = \frac{a\mathrm{g}_e \cos^2\phi + b\mathrm{g}_p\sin^2\phi}{\sqrt{a^2cos^2\phi + b^2\sin^2\phi}}\]

For numerical computations, a more convenient form is:

\[\mathrm{g} = \mathrm{g}_e\frac{1+k\sin^2\phi}{\sqrt{1-e^2\sin^2\phi}}\]

with the helper constant \(k\):

\[k = \frac{b\mathrm{g}_p}{a\mathrm{g}_e}-1\]

Parameters:

• lat (float) – Geographical latitude, in decimal degrees.

• h (float, default: 0.0) – Mean sea level height, in meters.

Returns:

g – Normal gravity at given point in space, in m/s^2.

Return type:

float

Examples

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

property normal_gravity_constant: float

Normal Gravity Formula Constant \(m\), computed as:

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

Example

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

property normal_gravity_potential: float

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

property polar_normal_gravity: float

Normal Gravity \(\mathrm{g}_p\) at the Pole, in \(\frac{\mathrm{m}}{\mathrm{s}^2}\), computed as:

\[\mathrm{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

property second_degree_zonal_harmonic: float

WGS 84 Second Degree Zonal Harmonic \(C_{2,0}\), computed as:

\[C_{2,0} = -\frac{J_2}{\sqrt{5}}\]

Example

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

property second_eccentricity_squared: float

Second Eccentricity Squared \(e'^2\), computed as:

\[e'^2 = \frac{a^2-b^2}{b^2}\]

Example

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

property sidereal_day: float

Sidereal Day, \(T_{sid}\), in seconds, computed as:

\[T_{sid} = \frac{2\pi}{\omega}\]

Example

>>> wgs = ahrs.utils.WGS()
>>> wgs.sidereal_day
86164.090530833

vertical_curvature_radius(lat: float) → float

Radius of the curvature in the prime vertical, estimated at a given latitude, \(\phi\), 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 the normal gravity, \(\mathrm{g}\), using 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, at a certain latitude, \(\phi\), can be written in the form of a series:

\[\mathrm{g} = \mathrm{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}\]

and \(\mathrm{g}_e\) is the measured normal gravity on the Equator. 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:

\[\mathrm{g} = \mathrm{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:

\[\mathrm{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\), \(\mathrm{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	\(\mathrm{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'.

Returns:

g – Normal gravity, in m/s^2, at given latitude.

Return type:

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 variations of gravity [WELMEC2021].

Manufacturers may adjust their instruments using the reference gravity formula:

\[\mathrm{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 above sea level in meters.

Parameters:

• lat (float) – Geographical Latitude, in decimal degrees.

• h (float, default: 0.0) – Mean sea level height, in meters.

Returns:

g – Normal gravity at given point in space, in m/s^2.

Return type:

float

Examples

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