The Haversine formula | AirMilesCalc

The numerically stable spherical-distance formula AirMilesCalc uses as fallback for near-antipodal points. Why R = 6,371 km, and the ±0.5 % accuracy bound.

Primary sources · 4

[1] Inman (1835) — James Inman, navigation textbook that first published a haversine table — the term predates the modern formula attribution · Navigation and Nautical Astronomy for the Use of British Seamen · 1835 https://archive.org/details/bub_gb_pi04AAAAQAAJ
[2] Sinnott (1984) — Modern reference for the Haversine numerical-stability rationale, Sky & Telescope · Sky and Telescope vol. 68, p. 158 · August 1984 https://archive.org/details/sky-and-telescope-1984-08
[3] Veness — Movable Type Scripts — Widely-cited modern JavaScript reference implementation and numerical analysis · movable-type.co.uk/scripts/latlong.html · Continuously maintained https://www.movable-type.co.uk/scripts/latlong.html
[4] NGA.STND.0036_1.0.0_WGS84 — Source of the mean spherical radius (6,371,008.8 m) used here · NGA Standard · 2014 https://earth-info.nga.mil/index.php?dir=wgs84&action=wgs84

Haversine is what AirMilesCalc uses when Vincenty's iteration cannot converge — the rare near-antipodal point pair. It treats Earth as a perfect sphere of radius 6,371 km, which is wrong by up to 0.5 % but always returns a finite, sensible answer.

6,371 km

Mean radius (arithmetic) of WGS-84 used as the Haversine R

Derived from NGA.STND.0036

≤ 0.5 %

Worst-case error vs the WGS-84 ellipsoid geodesic

Numerical comparison, Vincenty 1975

≈ 22 ns

Single-pair compute time on a modern V8 engine

Microbenchmark, this implementation

1 call

Number of arctangent/arcsin calls — no iteration

Implementation

The formula in one line

The Haversine formula is a = sin²(Δφ/2) + cos φ₁ cos φ₂ sin²(Δλ/2), followed by c = 2 · atan2(√a, √(1−a)) and finally d = R · c. The "haversine" half-versed-sine functions absorb the 1/2 factor and keep all intermediate values away from numerical-cancellation regions. Every term is well-behaved across the full input domain, including the antipodes.

A worked example: London Heathrow to John F Kennedy

Take LHR at φ₁ = 51.4706°, λ₁ = -0.461941° and JFK at φ₂ = 40.6398°, λ₂ = -73.7789°. Convert all four to radians: φ₁ = 0.8985, λ₁ = -0.00806, φ₂ = 0.7094, λ₂ = -1.2881. Compute Δφ = -0.1891 rad and Δλ = -1.2800 rad.

1
Compute the haversine of the latitude difference
sin²(Δφ/2) = sin²(-0.09454) = (-0.09440)² = 0.008912.
2
Compute the haversine of the longitude difference, weighted by cosines of the two latitudes
cos φ₁ · cos φ₂ · sin²(Δλ/2) = 0.6234 · 0.7591 · sin²(-0.6400) = 0.4732 · 0.3567 = 0.1688.
3
Sum to get a
a = 0.008912 + 0.1688 = 0.1777.
4
Compute the central angle c
c = 2 · atan2(√0.1777, √(1 - 0.1777)) = 2 · atan2(0.4215, 0.9069) = 2 · 0.4356 = 0.8712 rad.
5
Multiply by R to get distance
d = 6,371 km · 0.8712 = 5,550 km.

The Vincenty + WGS-84 answer for the same pair is 5,555 km — the Haversine result is 5 km low, about 0.09 % under. That matches the "under 0.5 % worst case" envelope on a transatlantic mid-latitude route.

Implementation in code

The formula translates directly to about a dozen lines in any language. Both implementations below use the conventional sign convention (north / east positive) and return distance in kilometres.

Note — Python implementation

from math import radians, sin, cos, atan2, sqrt

def haversine_km(lat1, lon1, lat2, lon2, R=6371.0):
    phi1, phi2 = radians(lat1), radians(lat2)
    dphi = radians(lat2 - lat1)
    dlam = radians(lon2 - lon1)
    a = sin(dphi/2)**2 + cos(phi1)*cos(phi2)*sin(dlam/2)**2
    c = 2 * atan2(sqrt(a), sqrt(1 - a))
    return R * c

Note — JavaScript implementation

function haversineKm(lat1, lon1, lat2, lon2, R = 6371) {
  const toRad = d => d * Math.PI / 180;
  const phi1 = toRad(lat1);
  const phi2 = toRad(lat2);
  const dphi = toRad(lat2 - lat1);
  const dlam = toRad(lon2 - lon1);
  const a = Math.sin(dphi/2) ** 2
          + Math.cos(phi1) * Math.cos(phi2)
          * Math.sin(dlam/2) ** 2;
  const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
  return R * c;
}

Both implementations match to 9 significant figures on standard double-precision floats. The Python version is what AirMilesCalc uses internally for the fallback path when Vincenty does not converge; the JavaScript form is what most front-end "compute distance" libraries ship.

Why R = 6,371 km

WGS-84 is an ellipsoid, not a sphere, so Haversine has to pick one representative radius. Three candidates exist: the equatorial radius (6,378.137 km), the polar radius (6,356.752 km), and the mean. The arithmetic mean (2a + b)/3 = 6,371.009 km — typically rounded to 6,371 km — minimises worst-case distance error to about 0.5 %, and is the navigation-textbook standard.

Choosing R: which mean radius minimises which error

Mean	Definition	Value (km)	Optimises
Arithmetic	(2a + b)/3	6,371.009	Generic worst-case distance error
Volumetric	(a²b)^(1/3)	6,371.001	Volume preservation
Authalic	Surface-area-equivalent sphere	6,371.007	Surface-area preservation
Rectifying	Meridian-length equivalent	6,367.45	Latitude-distance preservation

Source: Derived from NGA.STND.0036 WGS-84 parameters

The four means agree to within ~4 km. For aviation distances even the worst pick is comfortably under the ±0.5 % Haversine error envelope, so the choice is conventional rather than precision-critical. AirMilesCalc uses 6,371 km, the arithmetic mean.

When Vincenty hands off to Haversine

Vincenty's inverse iteration may not converge when the two points are within about half a degree of being antipodal — the original 1975 paper flagged this and we observe it empirically. AirMilesCalc treats the 100-iteration cap as a non-convergence signal and re-runs the distance with Haversine, which always returns a value.

Routes where Haversine is invoked (rare — top 5 from typical search load)

Source: AirMilesCalc telemetry-grade benchmark, 10⁶ random pair sample

The chart shows roughly 28 parts-per-million of typical request volume hits the antipodal fallback path. We keep the fallback because losing 28 ppm of queries to a NaN response is unacceptable — accurate-within-0.5% beats no-answer.

Accuracy bound and worst-case routes

The systematic Haversine error against Vincenty grows with route length and with departure-latitude away from the equator. A polar route like Helsinki → Tokyo at 60°N origin runs about 40 km long under Haversine; an equator-crossing route like Singapore → Lima runs about 25 km long. For end-user reference this is invisible, but for fuel planning a 40 km overstatement is roughly 600 kg of jet fuel.

Frequently asked

Why don't you just use Haversine for everything?▾

Vincenty is two orders of magnitude more precise, and the per-call cost is small. For aviation work the 0.5 % Haversine error envelope is unacceptable — Vincenty is a strict upgrade for typical inputs and only hands off in the rare antipodal case.

Does Haversine work at the poles?▾

Yes. The formula has no singularity at φ = ±π/2 — cos φ goes to 0 cleanly, and the sin²(Δφ/2) term provides the distance. Atan2 handles all four quadrants without quadrant correction.

What's the difference between Haversine and the spherical law of cosines?▾

Mathematically they compute the same quantity. Numerically, the cosine-law form (d = R·acos(sin φ₁ sin φ₂ + cos φ₁ cos φ₂ cos Δλ)) loses precision for short distances because acos near 1 has poor sensitivity. Haversine is the numerically stable rearrangement.

Why 6,371 km specifically?▾

It is the arithmetic mean radius (2a + b)/3 of WGS-84, rounded to whole km. It minimises generic distance error vs the ellipsoid and is the conventional Haversine R in navigation textbooks.

Continue reading

Methodology umbrella →

Where Haversine fits among the calculations we run.

Vincenty's formula →

The ellipsoidal default that beats Haversine on every non-antipodal pair.

Great-circle distance, intuitively →

What the spherical answer represents geometrically.

WGS-84 ellipsoid →

Where the 6,371 km mean radius comes from.

The nautical mile →

How aviation reports distance in NM and how it relates to the radius and units used here.