WGS84 Gravity Model — Robotics/Gravity Models

WGS84 Gravity Model

Robotics / Gravity Models

Normal gravity of the 1984 World Geodetic System's equipotential ellipsoid (NIMA TR8350.2) at a geodetic position, in the local north-east-down frame: g = [g_N; 0; g_D], with g_D positive downward and g_N positive northward. With the geodetic latitude φ, the height h and the ellipsoid's constants a, e², f, γ_e, k:

• Taylor Series – g_D = γ(φ)·(1 − 2(1 + f + m − 2f·sin²φ)·h/a + 3h²/a²), with γ(φ) = γ_e(1 + k·sin²φ) / √(1 − e²sin²φ); g_N = 0.
• Close Approximation – the point is taken to ellipsoidal-harmonic coordinates (u, β), the two components γ_u, γ_β of the normal field are evaluated there, and g_D = √(γ_u² + γ_β²); g_N = 0.
• Exact – the same two components rotated through the geocentric latitude ψ and the deflection α = φ − ψ into the geodetic vertical (g_D) and the north tangent (g_N), so the north component is real.

Ports

• lla – the position as one [3,1], [φ; λ; h]: geodetic latitude and longitude in degrees, then height above the ellipsoid in metres, in that order. A latitude past a pole is folded back and its longitude turned half a revolution, as the reference does.
• g – the gravity vector as one [3,1], [g_N; 0; g_D] in m/s². The east entry is always 0, and so is the north entry except under Exact.

Parameters

• Gravity Model – which method runs. This selects which arithmetic executes, not just a value:
  • Taylor Series – the closed-form surface gravity expanded to second order in h. The default, and Simulink's. Meant for low heights; the three switches below do not apply to it.
  • Close Approximation – the ellipsoidal-harmonic field at the point, as a magnitude. Sub-microgal below about 20000 m.
  • Exact – the same field resolved into down and north.
• Exclude Atmosphere – on (the default) uses GM without the atmosphere's mass, 3.9860009×10¹⁴; off uses the WGS84 value including it, 3.986004418×10¹⁴ m³/s². Harmonic methods only.
• Precessing Reference Frame – on (the default) takes Earth's rate as the IAU value plus the precession in right ascension at the date below, ω = 7.2921151467×10⁻⁵ + 7.086×10⁻¹² + 4.3×10⁻¹⁵·T, T in Julian centuries from J2000; off takes the WGS84 rate 7.292115×10⁻⁵ rad/s. Harmonic methods only.
• No Centrifugal Effects – on (the default) removes the centrifugal term, leaving pure attraction; off folds it in, which is what an instrument at rest on the Earth feels. Harmonic methods only.
• Month – the date's month, January to December. Defaults to October. Read only with a precessing frame.
• Day – the date's day, a whole number from 1 to 31. Defaults to 10. A day past the month's end is taken as its last day (February 30 as February 28), as Simulink does.
• Year – the date's year, a whole number of 1 or more. Defaults to 2004. A year of 99 is special, exactly as in the reference: the Julian date is then the Day itself rather than a date in the year 99.
• Sampling Time (s) – zero or less inherits the solver's rate; a positive value runs the block at that period.

Code export

All ten targets: Python, MATLAB, Java, Rust, C, C++, VHDL, Verilog, SystemVerilog and PLC Structured Text.

The method and the three switches are baked in at export time: a core built for Taylor Series contains no harmonic arithmetic at all, and one built without the centrifugal term contains no centrifugal term. Earth's rate is computed from the date once, when the block's configuration is loaded, and emitted as a constant. Every target emits the reference's association operation for operation, because two of the intermediate quantities are small differences of large ones and a reordered product moves the last three digits; the C and C++ bodies also switch off fused multiply-add for the same reason.

The three hardware targets are simulation-only: sines, cosines, arctangents, square roots and divisions of a planet-sized geometry do not belong in a Q16.16 datapath, so the arithmetic runs in real and only the ports are fixed point. Q16.16 also bounds what a port can carry to about ±32767, which keeps a hardware core below about 32 km of height. VHDL carries the body as a function with its own real variables; PLC Structured Text has no floor and rebuilds it from TRUNC.

How closely a target agrees depends on its math library, because two of the intermediate quantities are small differences of large ones and magnify a rounding about a hundred thousand times. Measured over 16 configurations: the generated MATLAB reproduces the Simulink block exactly, and C, Python, Verilog and SystemVerilog (under Icarus) agree to about 5×10⁻¹⁵ m/s². A library that differs in the last place moves the answer by up to about 2×10⁻¹²: Java's does (7×10⁻¹³), and so does an optimising C or C++ compiler that fuses a sine and a cosine of the same angle into one call. VHDL under GHDL uses the IEEE reference MATH_REAL, whose sine and arctangent are good to about 10⁻⁸, and there agrees to about 5×10⁻⁴ m/s²; its cosine of π/2 is exactly zero, which would divide by zero at a latitude of exactly ±90°, so the VHDL body replaces that zero with the correctly rounded value.

Simulink bridge

Import and export, mapped to Aerospace Blockset's aerolibgravity2/WGS84 Gravity Model   – the library name really ends in two spaces. Gravity Model → model (Taylor Series/Close Approximation/Exact → WGS84 Taylor Series/WGS84 Close Approximation/WGS84 Exact), Exclude Atmosphere → no_atmos, Precessing Reference Frame → precessing, No Centrifugal Effects → no_centrifugal (all three on/off 1:1), Month → month (1:1), Day → day and Year → year. Every mapping is 1:1 and lossless in both directions.

Three Simulink parameters are always implied and carry no configuration here: units is always Metric (MKS) (English units do not cross), jd_loc is always off – Simulink's option of taking the Julian date on a second input port does not cross, because no configuration here can add a port, and a Julian date near 2.46×10⁶ could not ride a hardware port in any case – and action is always Warning. That last one only chooses whether Simulink reports a height below 0 or above 20000 m, or a folded latitude; the values are the same either way, and this block computes them without a report. An imported block holding any other value of the three is reported rather than silently accepted. Simulink refuses a [3,1] MATRIX on this port: it takes a 1-D 3-element vector (or an [m,3] array of points), so feed the exported block from something that produces a vector. The Simulink block defines no SampleTime, so the rate stays on the ICore side.

Notes

• Algebraic and stateless: the output depends only on this sample's input.
• Not linear, so the block carries no state space and model reduction correctly reports it as unmergeable.
• The longitude does not change the answer – the ellipsoid is symmetric about the spin axis – except in its last digits, where the reference's own rounding of it shows. It is wrapped and folded exactly as the reference does so that those digits agree too.
• A latitude past a pole is folded in radians (π − φ across the north pole, −π − φ across the south), with the longitude turned by +π and −π respectively. Only Exact can tell the difference, and without the fold a latitude of 181 is wrong by 19.6 m/s² there.
• Do not check Close Approximation or Exact against MATLAB's gravitywgs84 at its defaults: the function includes the centrifugal term and the atmosphere by default and the Simulink block (and this one) excludes both, so the two disagree by about 0.017 m/s² at 45°.
• Verified against R2026a: over 16 configurations and 25 points (both hemispheres, latitudes of 95, 181, 270 and 450, heights from −500 to 32000 m, the leap-day, invalid-day and year-99 dates) the Simulink block and this arithmetic agree to 5.3×10⁻¹⁵.