Module geodesic 

Geodesics on an ellipsoid: the distance and bearing between two places, and the place a given distance away along a bearing.

A geodesic is the shortest path between two points on the ellipsoid’s surface. The two problems are Karney’s:

• inverse (Ellipsoid::geodesic_inverse): given two points, the length s₁₂ of the geodesic between them and its azimuths α₁ (leaving point 1) and α₂ (arriving at point 2);
• direct (Ellipsoid::geodesic_direct): given a point, an azimuth α₁ and a distance s₁₂, the end point and α₂.

The algorithms are C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87 (2013) 43–55, doi:10.1007/s00190-012-0578-z (arXiv:1109.4448v2). The ellipsoid is mapped onto an auxiliary sphere, where a geodesic is a great circle, and distance and longitude are corrected by series in the third flattening n = f/(2 − f) to sixth order in the flattening f (§2, the direct problem); the inverse finds α₁ by Newton’s method (§4), from a starting guess (§5). The code is GeographicLib’s, as georust’s geographiclib-rs ports it (MIT). Karney states that round-off in both problems stays under 15 nanometers on WGS 84 (§7, page 10), and that for f up to 1/150 the series’ truncation is smaller than round-off (page 9): past that the series lose accuracy, so a flatter ellipsoid is refused (GEODESIC_MAX_FLATTENING). Karney’s published test set of 500,000 WGS 84 geodesics (doi:10.5281/zenodo.32156) is the check: tests/geodtest.rs, with the measured errors in docs/physics/geodesy.md.

Azimuths are clockwise from north, in radians, in [−π, π]; .rem_euclid(TAU) gives the 0 to 2π of a compass. α₂ is the direction of travel at point 2, so the bearing back to point 1 from there is α₂ ± π. Heights play no part: the path lies on the ellipsoid’s surface, not at the points’ heights.

Structs

GeodesicDirect

The end of a geodesic from a start, an azimuth and a distance: Karney’s direct problem.

GeodesicInverse

The geodesic between two points: Karney’s inverse problem.

Constants

GEODESIC_MAX_FLATTENING

The largest flattening the geodesics accept, 1/150: Karney 2013 (page 9) shows the sixth-order series’ truncation below f64 round-off up to there. GeographicLib’s own error table for the series grows to 10 µm at f = 0.05 and 0.3 m at 0.2. Every planet-like body hpr flies on is inside: WGS 84’s f is 1/298.257.