geo.distance

export function distance(fromLat: number, fromLng: number, toLat: number, toLng: number): number

import { distance } from "#fune/geo.distance@^2";

import { atan } from "./math_atan.ts";  ← from math.atan ^1.0.0 · built alongside by fune
import { roundFloat } from "./math_round_float.ts";  ← from math.round-float ^1.0.0 · built alongside by fune
import { sinCos } from "./math_sin_cos.ts";  ← from math.sin-cos ^1.0.0 · built alongside by fune

/**
 * Great-circle distance between two points, by the haversine formula, on a
 * sphere of the IUGG mean Earth radius, in metres rounded to millimetres.
 *
 * Sine, cosine and arctangent come from math.sin-cos and math.atan, not
 * Math.sin or Math.atan2: those come from each platform's maths library, which
 * may differ in the last bit between JavaScript, Python and Rust, and a value
 * next to a rounding boundary can then go either way. With the shared helpers
 * the unrounded distance is the same double in every language, and
 * math.round-float rounds that double the same way everywhere.
 */

// IUGG mean radius R1 = (2a + b) / 3 of the GRS 80 ellipsoid (Moritz,
// "Geodetic Reference System 1980", Journal of Geodesy 74 (2000) 128-133).
const EARTH_RADIUS_METRES = 6371008.8;

const PI = 3.141592653589793;
const HALF_PI = PI / 2;
const DEGREES = PI / 180;

// atan2 for y >= 0 and x >= 0, which is all haversine needs. x is 0 only when
// the points are exactly antipodal (a = 1); y / 0 would be infinite, which
// math.atan refuses.
function atan2Positive(y: number, x: number): number {
  if (x === 0) return y === 0 ? 0 : HALF_PI;
  return atan(y / x);
}

function checkLatitude(value: number): void {
  if (typeof value !== "number" || !Number.isFinite(value)) {
    throw new TypeError("latitude must be a finite number of degrees");
  }
  if (value < -90 || value > 90) throw new RangeError("latitude must be between -90 and 90 degrees");
}

function checkLongitude(value: number): void {
  if (typeof value !== "number" || !Number.isFinite(value)) {
    throw new TypeError("longitude must be a finite number of degrees");
  }
  if (value < -180 || value > 180) throw new RangeError("longitude must be between -180 and 180 degrees");
}

export function distance(fromLat: number, fromLng: number, toLat: number, toLng: number): number {
  checkLatitude(fromLat);
  checkLongitude(fromLng);
  checkLatitude(toLat);
  checkLongitude(toLng);

  // A longitude difference of 359 degrees is a 1 degree step across the
  // antimeridian; sin^2 of the half-angle is the same either way, so no
  // wrapping is needed.
  const sinHalfLat = sinCos(((toLat - fromLat) * DEGREES) / 2).sin;
  const sinHalfLng = sinCos(((toLng - fromLng) * DEGREES) / 2).sin;
  const cos1 = sinCos(fromLat * DEGREES).cos;
  const cos2 = sinCos(toLat * DEGREES).cos;

  let a = sinHalfLat * sinHalfLat + cos1 * cos2 * sinHalfLng * sinHalfLng;
  // Rounding can push a a hair outside [0, 1] for antipodal points.
  if (a < 0) a = 0;
  if (a > 1) a = 1;
  // atan2 rather than asin(sqrt(a)): stays accurate for antipodal points too.
  const c = 2 * atan2Positive(Math.sqrt(a), Math.sqrt(1 - a));
  // Half away from zero on the exact value of the double (math.round-float).
  return roundFloat(EARTH_RADIUS_METRES * c, 3);
}

Install

fune build

With that line in your source, in a TypeScript project (language typescript in fune.project), fune build resolves it and its 3 dependencies, pins them in fune.lock, downloads only the TypeScript package of each, and builds the code above into your project’s .fune/build, one readable file per capability with a header linking back here. Or pin a range in fune.project and build in one step:

fune add geo.distance

The manifest, vectors and README with only the TypeScript implementation. Install it without the registry with fune add ./geo.distance-2.0.0-typescript.fune, or fetch it from a terminal with fune pull geo.distance@2.0.0:typescript.

The whole function, every language, is one file too: geo.distance-2.0.0.fune, 20,447 bytes, sha256 2b5c24b572d5c3696eb23324dcd7c6dc2fd9d26a25655445ec770267aa124a4d. It installs into a project of any language.

// fune: before geo.distance

// fune: after geo.distance

// fune: replace math.atan in geo.distance
// fune: replace math.round-float in geo.distance
// fune: replace math.sin-cos in geo.distance

// fune: step geo.distance after <n|label>

More from the author

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.

Files