# geo.distance

The great-circle distance between two latitude/longitude points, in metres,
by the haversine formula, rounded to the millimetre.

## The model

The Earth is taken as a sphere of radius 6,371,008.8 m, the IUGG mean radius
R1 = (2a + b) / 3 of the GRS 80 ellipsoid (H. Moritz, "Geodetic Reference
System 1980", Journal of Geodesy 74 (2000) 128-133; also Bulletin Géodésique
54 (1980)). Many snippets use 6,371,000 m or the equatorial 6,378,137 m; those
give different answers, which is why the radius is part of the contract.

A sphere is not the Earth. Against the WGS 84 ellipsoid (Vincenty, Karney) the
haversine distance is off by up to about 0.5%. Use this for "how far is the
nearest depot", delivery radii and sorting by distance, not for surveying.
The millimetre rounding is about reproducibility, not accuracy.

Longitudes do not need wrapping: 179.5 to -179.5 is 1 degree across the
antimeridian, because the formula only uses sin² of half the difference.
Coordinates outside -90..90 and -180..180 are errors rather than being wrapped,
since they usually mean latitude and longitude were swapped.

The result uses atan2(√a, √(1−a)) rather than asin(√a). For points closer than
a metre or so it keeps millimetre accuracy, where the spherical law of cosines
(acos) loses it. Near-antipodal points are ill-conditioned for any haversine in
double precision: within about a degree of the antipode the answer can be a
few millimetres from the exact spherical value (elsewhere it is within 1e-7 m).

## Why the answers are identical in every language

JavaScript's Math.sin, Python's math.sin and Rust's f64::sin come from
different maths libraries, which are allowed to differ in the last bit.
Rounding to millimetres hides that almost always, but not always: a value that
lands on a rounding boundary can go either way.

So this capability does not call them. It carries its own sin, cos and atan,
built only from operations IEEE 754 requires to be correctly rounded
(+ − × ÷, sqrt, floor), evaluated in the same order in all three languages:
range reduction to a quarter turn, fixed-length Taylor polynomials in Horner
form, and angle halving for atan. No language fuses a multiply and an add
behind your back (JavaScript and Python cannot; Rust does not without an
explicit mul_add). The unrounded distance is therefore the same double in
TypeScript, Python and Rust, and the rounding step is only presentation.
This was checked on 4,000 random pairs (including near-coincident and
near-antipodal ones): the unrounded and rounded results were bit-for-bit
identical in all three.

## Rounding

Half away from zero, applied to the binary64 value: y = d × 1000,
r = floor(y), r + 1 if y − r ≥ 0.5, result r / 1000. Not Python's round(),
which rounds halves to even.