# 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 (see Accuracy).
## What changed in 2.0.0
- **Rounding is on the exact value of the double.** 1.x multiplied by 1000
and rounded that product, so a distance stored just below a millimetre tie
could round up: the double 14704080.96649999916... (Tokyo Haneda,
35.5494, 139.7798, to -33.9399, 18.786034 near Stellenbosch) times 1000 rounds to exactly
14704080966.5, and 1.x said 14704080.967. 2.0.0 rounds with
`math.round-float`, half away from zero on the value actually held, and says
14704080.966, the way 2.675 to two places is now 2.67 rather than 2.68. The
high-precision reference agrees: the true distance for those inputs is
14704080.96649999826... m, below the tie as well.
- **Sine, cosine and arctangent come from `math.sin-cos` and `math.atan`**
instead of private copies. `math.sin-cos` reduces the angle with a two-part
π/2 and `math.atan` is fdlibm's arctangent, both from + − × ÷ only, so the
unrounded distance can differ from 1.x in the last bits (by up to about
1e-9 m on the vectors) while staying the same double in every language.
The one piece left here is a four-line atan2 for y, x ≥ 0 on top of
`math.atan` (atan(y/x), and π/2 when x is 0, which only exactly antipodal
points reach), since the registry has no atan2.
- **Vectors.** Every 1.x vector keeps its expected answer; each was re-derived
from an independent 60-digit reference (own series for π, sin, cos and
atan, decimal square root), and the rounded true distance equals it. Two
vectors were added that show the rounding change: Tokyo Haneda to near
Stellenbosch (1.x .967, now .966) and London to 33.8688°S 151.714482°E in
the Tasman Sea (17024346.35449999943..., 1.x .355, now .354). The vectors
now compare exactly (`floats exact`): the three languages return the same
double.
Callers of 1.x who want the new rounding should move to ^2; 1.x is unchanged.
## 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. Its trigonometry comes from
`math.sin-cos` and `math.atan`, built only from operations IEEE 754 requires
to be correctly rounded, in the same order in every language; the square roots
here are IEEE's correctly rounded sqrt. 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 `math.round-float` rounds it the same way in
all three.
## Accuracy
Measured against the 60-digit reference, with the inputs taken as the exact
doubles given, on 5,500 uniformly random pairs plus 1,500 of each special kind
(coordinates to 6 decimals):
- anywhere on the globe: the unrounded distance is within 4e-8 m of the true
spherical distance (the largest errors are on the longest distances; at
20,000 km that is about ten units in the last place);
- points less than about 150 m apart: within 5e-14 m;
- 1° to 7° from the antipode: within 3e-7 m; within 1° of the antipode:
within 3e-6 m, because the haversine itself is ill-conditioned there.
1.x had the same error bounds; they come from converting degrees to radians
and from the formula's conditioning, not from the sine or arctangent. On every
sampled pair the rounded answer equalled the true distance rounded half away
from zero (1.x missed one, near the antipode). It cannot be guaranteed: when
the true distance is within a few nanometres of a millimetre tie, the double
can sit on the other side of it (35.5494, 139.7798 to -33.9399, 18.933865 is
truly 14691235.9925000015... m, but its double is 14691235.99249999970... m,
so the answer is 14691235.992). The rounding is exact on the double; the
double is within the bounds above of the truth.
## Rounding
Half away from zero to 3 decimals, by `math.round-float`, on the exact value
of the double: the nearest multiple of 0.001 to the value held, an exact tie
going away from zero. Not Python's round(), which rounds halves to even, and
not floor(d × 1000 + 0.5), which rounds the product first. A result is never
-0.