Geodesic distance (haversine, Vincenty, Karney)

Deep dive

How it works

In: two (lat, lon) points in degrees, an ellipsoid (WGS84) or a sphere radius · Out: distance in metres (and azimuths for Karney) · State: none · Cost: O(1) per pair; order of 28 ns haversine, 237 ns Vincenty (measured below)

1. Sphere (haversine). h = sin²(Δφ/2) + cos φ1 cos φ2 sin²(Δλ/2), d = 2R asin(√h). The half-angle form stays accurate for tiny distances, where the law of cosines (acos) loses digits.
2. Chord shortcut. Map each point to a unit vector once. Euclidean chord length is monotone in great-circle distance, so sort, threshold and k-d-tree search on chords and convert only the survivors: d = 2R asin(chord/2).
3. Ellipsoid (Vincenty). Use reduced latitudes U = atan((1-f) tan φ). Iterate the longitude λ on the auxiliary sphere until it moves less than 1e-12 rad (about 0.06 mm), then apply the series in u² = cos²α (a²-b²)/b² to turn the auxiliary arc σ into metres.
4. Ellipsoid (Karney). Same idea, but the inverse problem is solved by Newton's method with bisection as a fallback, so every input converges. Use a library.

Use WGS84: a = 6378137 m, f = 1/298.257223563. The mean radius (2a+b)/3 is 6371008.8 m. The sphere ignores the 21 km difference between equatorial and polar radius.

Code

use std::f64::consts::PI;

const A: f64 = 6_378_137.0;            // WGS84 semi-major axis, m
const F: f64 = 1.0 / 298.257_223_563;  // WGS84 flattening
const B: f64 = A * (1.0 - F);
const R_MEAN: f64 = 6_371_008.8;       // (2a+b)/3

/// Sphere, haversine form. Degrees in, metres out.
pub fn haversine(lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> f64 {
    let (p1, p2) = (lat1.to_radians(), lat2.to_radians());
    let (dp, dl) = ((lat2 - lat1).to_radians(), (lon2 - lon1).to_radians());
    let h = (dp / 2.0).sin().powi(2) + p1.cos() * p2.cos() * (dl / 2.0).sin().powi(2);
    2.0 * R_MEAN * h.sqrt().min(1.0).asin()
}

/// Unit vector; chord distance ranks like great-circle distance.
pub fn unit(lat: f64, lon: f64) -> [f64; 3] {
    let (p, l) = (lat.to_radians(), lon.to_radians());
    [p.cos() * l.cos(), p.cos() * l.sin(), p.sin()]
}
pub fn chord(a: [f64; 3], b: [f64; 3]) -> f64 {
    ((a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2) + (a[2] - b[2]).powi(2)).sqrt()
}
pub fn chord_to_m(c: f64) -> f64 { 2.0 * R_MEAN * (c / 2.0).min(1.0).asin() }

/// WGS84, Vincenty inverse. None if it fails to converge (near-antipodal).
pub fn vincenty(lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> Option<f64> {
    let l = ((lon2 - lon1 + 540.0).rem_euclid(360.0) - 180.0).to_radians(); // wrap to [-180,180)
    let u1 = ((1.0 - F) * lat1.to_radians().tan()).atan();
    let u2 = ((1.0 - F) * lat2.to_radians().tan()).atan();
    let (su1, cu1, su2, cu2) = (u1.sin(), u1.cos(), u2.sin(), u2.cos());
    let mut lam = l;
    let (mut ss, mut cs, mut sg, mut c2a, mut c2sm) = (0.0, 0.0, 0.0, 0.0, 0.0);
    for _ in 0..200 {
        let (sl, cl) = (lam.sin(), lam.cos());
        ss = ((cu2 * sl).powi(2) + (cu1 * su2 - su1 * cu2 * cl).powi(2)).sqrt();
        if ss == 0.0 { return Some(0.0); } // coincident points
        cs = su1 * su2 + cu1 * cu2 * cl;
        sg = ss.atan2(cs);
        let sa = cu1 * cu2 * sl / ss;
        c2a = 1.0 - sa * sa;
        c2sm = if c2a == 0.0 { 0.0 } else { cs - 2.0 * su1 * su2 / c2a }; // equatorial line
        let c = F / 16.0 * c2a * (4.0 + F * (4.0 - 3.0 * c2a));
        let next = l + (1.0 - c) * F * sa
            * (sg + c * ss * (c2sm + c * cs * (-1.0 + 2.0 * c2sm * c2sm)));
        let done = (next - lam).abs() < 1e-12; // about 0.06 mm
        lam = next;
        if done { break; }
        if lam.abs() > PI { return None; }
    }
    let u2s = c2a * (A * A - B * B) / (B * B);
    let ka = 1.0 + u2s / 16384.0 * (4096.0 + u2s * (-768.0 + u2s * (320.0 - 175.0 * u2s)));
    let kb = u2s / 1024.0 * (256.0 + u2s * (-128.0 + u2s * (74.0 - 47.0 * u2s)));
    let dsg = kb * ss * (c2sm + kb / 4.0 * (cs * (-1.0 + 2.0 * c2sm * c2sm)
        - kb / 6.0 * c2sm * (-3.0 + 4.0 * ss * ss) * (-3.0 + 4.0 * c2sm * c2sm)));
    Some(B * ka * (sg - dsg))
}

Checked: Flinders Peak to Buninyong (-37°57'03.72030" 144°25'29.52440" to -37°39'10.15610" 143°55'35.38390") gives 54972.271 m (Vincenty's test case; haversine 54925.5 m); one degree of equatorial longitude a·π/180 = 111319.4908 m; pole to equator 10001965.729 m; (0,0) to (0.5,179.5) = 19936288.579 m; (0,0) to (0.5,179.7) returns None; the chord route equals haversine to 1 mm.

Complexity and constants

• Haversine: 2 sin, 2 cos, sqrt, asin; order of 28 ns per pair measured (Apple M-series, release build, 200k random pairs, scalar f64).
• Vincenty: 5-10 iterations, each sin + cos + atan2; order of 237 ns measured, 8x haversine. Karney's library is of similar order and robust.
• Chord distance after the unit vectors exist: 3 subtractions, 3 multiplies, no transcendental function. Unit vectors cost 24 bytes per point as f64 (12 as f32).
• Error against Vincenty over random pairs: max 0.56%, also 0.56% for pairs within ~10 km. One degree of latitude is 110574 m at the equator and 111694 m at 89° (measured); R = 6371008.8 m has no way to represent that.
• Vincenty failed to converge for 2 of 200,000 random pairs (measured).

Efficiency

• Cost: one haversine on an Apple M-series core, scalar f64, order of 28 ns; one Vincenty inverse order of 237 ns (both measured here over 200,000 random pairs). Bound: compute (transcendentals); at millions of candidates, the cache misses of fetching the points dominate.
• Levers:
  • Prefilter before any trigonometry — a cell index range scan or a k-d tree over unit vectors cuts a billion-pair job to the few candidates — always for radius and nearest queries.
  • Rank on chord² or on h, take asin/sqrt only for the final answer — saves one transcendental per candidate — top-k and radius filters.
  • Precompute per point (cos φ, unit vector, radians) in struct-of-arrays — lets the compiler vectorise the loop — any repeated query over a fixed point set.
  • Cheapest model under the error budget — sphere is about 8x cheaper than ellipsoid here (measured) — sphere for filters, ellipsoid only for the shortlist or the displayed value.
  • Project once and use Euclidean distance — no trigonometry per pair — inside one city (see map projections).
• Naive → standard → state of the art:
  • Euclidean on raw degrees · cheapest, wrong by the factor cos φ in longitude.
  • Haversine on a sphere · order of 30 ns, error up to about 0.5%.
  • Karney (2013) on WGS84 · same order as Vincenty, under 15 nm, always converges.

Implementation notes

• Wrap the longitude difference; points on both sides of ±180° are close. The code wraps l in Vincenty; haversine is unaffected because sin²(Δλ/2) is periodic.
• Clamp h.sqrt() to 1 before asin: rounding can push it to 1.0000000000000002.
• Near-antipodal points make the haversine asin ill-conditioned (error grows when h approaches 1); use atan2(√h, √(1-h)) if that matters.
• Keep latitudes inside [-90, 90]; tan(90°) in Vincenty is large but the atan of it is fine in f64.
• Use f64: f32 has about 7 digits, so latitude alone resolves only order of 0.4 m and 10,000 km distances round to order of 1 m.
• If you need the path, azimuths or points along the path, call Karney's direct/inverse and GeodesicLine.

Pitfalls

• Using acos(sin sin + cos cos cos): catastrophic cancellation below ~1 km.
• Mixing radii: 6371 km (mean) against 6378.137 km (equatorial) is a 0.1% difference.
• Haversine for geofences at cm accuracy: the model error is hundreds of metres per 100 km.
• Vincenty on user input: antipodal failures turn into infinite loops if the iteration is unbounded.
• Degrees where radians are expected, and lat/lon argument order against lon/lat (GeoJSON, PostGIS) order.
• Great circles are not constant-bearing lines; rhumb line distance is a different, longer quantity.

Variants and alternatives

In production

• Databases and search engines: Redis GEODIST is spherical, its docs state errors up to 0.5%; Lucene's SloppyMath.haversinMeters is the haversine used for geo-distance sorting (docs).
• Mobile and navigation: Android Location.distanceBetween defines distance and bearing on the WGS84 ellipsoid and warns that near-antipodal points may give meaningless results (docs).
• Geodesy and surveying: the geod utility of PROJ (docs) and GeographicLib solve the direct and inverse problems.
• PostGIS ST_Distance on geography defaults to a spheroid distance (use_spheroid = true).
• PROJ's geod and GeographicLib implement Karney; geopy's geodesic and pyproj's Geod use it.

Implementations

Sources

• Charles F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43-55, 2013. paper · arXiv · code
• GeographicLib, Geodesics on an ellipsoid of revolution. docs
• T. Vincenty, Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations, Survey Review 23(176), 88-93, 1975.
• R. W. Sinnott, Virtues of the Haversine, Sky and Telescope 68(2), 159, 1984.
• Redis, GEODIST. docs · PostGIS, ST_Distance. docs