Primary sources · 4
- [1] Vincenty (1975a) — Direct and Inverse Solutions of Geodesics on the Ellipsoid with Application of Nested Equations · Survey Review XXIII (176), pp. 88–93 · April 1975 https://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
- [2] Vincenty (1975b, unpublished) — Alternative iteration scheme for near-antipodal points, referenced by Vincenty himself · Defense Mapping Agency Aerospace Center internal report · 1975 https://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
- [3] NGA.STND.0036_1.0.0_WGS84 — WGS-84 ellipsoid defining parameters (a, 1/f, GM, ω) · National Geospatial-Intelligence Agency Standard, v1.0.0 · July 2014 https://earth-info.nga.mil/index.php?dir=wgs84&action=wgs84
- [4] Karney (2013) — Algorithms for geodesics — a modern alternative to Vincenty with guaranteed convergence · Journal of Geodesy 87(1), pp. 43–55 · January 2013 https://doi.org/10.1007/s00190-012-0578-z
Vincenty's 1975 formula solves a problem that has no closed-form answer: what is the shortest path between two points on an oblate spheroid? It does so iteratively, converging on the answer through repeated trigonometric refinement, and lands within half a millimetre of the truth.
The problem Vincenty solved
A geodesic on a sphere is just an arc of a great circle, with a one-line formula. On an oblate spheroid the geodesic is a more complicated curve — it no longer lies in a single plane through the centre, and the latitude-as-a- function-of-longitude relationship has no elementary solution. Vincenty turned the problem into a fixed-point iteration on a single auxiliary longitude.
The five quantities the loop refines
The iteration moves a single number, λ, until successive estimates differ by less than the tolerance. λ is the longitude difference projected onto the auxiliary sphere of radius b/a, where b is the polar semi-axis and a the equatorial. Each pass computes σ (angular distance), sin α (azimuth-related), cos 2σm (the angular midpoint), and the correction term C.
| Symbol | Meaning | Range |
|---|---|---|
| σ (sigma) | Angular separation on the auxiliary sphere | 0 … π rad |
| sin α | Sine of the azimuth at the equator crossing of the geodesic | −1 … +1 |
| cos 2σm | Cosine of twice the angular distance from the equator crossing to the midpoint | −1 … +1 |
| C | Correction coefficient that captures ellipsoidal flattening | small, ≈ f scale |
| λ (lambda) | Longitude on the auxiliary sphere — the iteration's fixed point | −π … +π rad |
The convergence behaviour is quadratic for well-conditioned inputs. After two or three passes the residual is already below 10⁻⁴ rad and each subsequent pass squares the remaining error. For a typical commercial routing pair the loop exits in five or six passes — the 100-iteration cap protects only the pathological cases.
The four-step algorithm
- 1Reduce the latitudes
Convert each geographic latitude φ to a reduced latitude U via tan U = (1 − f) tan φ. This maps the ellipsoid problem onto a unit auxiliary sphere whose meridional radius depends on f. Cache sin U₁, cos U₁, sin U₂, cos U₂ — they appear in every subsequent expression.
- 2Iterate the longitude correction
Initialise λ to the geographic longitude difference L. Each pass computes sin σ, cos σ, then σ via atan2; then sin α and cos²α; then cos 2σm; then C. Finally update λ from the right-hand side of equation (14). Stop when |λ − λ_prev| < 10⁻¹² rad or after 100 passes.
- 3Detect non-convergence
Vincenty himself noted in the 1975 paper that "the formula will not converge in the case of nearly antipodal points." We treat the 100- iteration cap as a non-convergence signal and fall back to the Haversine spherical formula, which always returns a finite distance.
- 4Apply the ellipsoidal series correction
Compute u² = cos²α · (a² − b²) ⁄ b², then the series coefficients A and B (Vincenty equations (3), (4)). The surface distance is s = b · A · (σ − Δσ) where Δσ is the bracket of equation (6). Convert to km, statute miles (× 0.621371), or nautical miles (× 0.539957).
The antipodal edge case
When the two points are within roughly half a degree of being on opposite sides of the Earth, the iteration's fixed-point map loses its contraction property — λ can oscillate or even diverge. Vincenty published a follow-up report the same year with an alternative iteration that handles this case; because the report was internal and never widely typeset, most production implementations (ours included) use a Haversine fallback instead.
| Route | Distance (km) | Iterations | Behaviour |
|---|---|---|---|
| JFK → LHR | 5,555 | 5 | Smooth, residual drops 6 orders / pass |
| LHR → SYD | 17,016 | 6 | Smooth, large σ slows convergence one pass |
| Quito (UIO) → Padang (PDG, near-antipode) | ≈ 19,930 | fails @ 100 | Haversine fallback engaged |
| Same airport (φ₁=φ₂, λ₁=λ₂) | 0 | 0 | Short-circuit: sin σ = 0 → return 0 |
Comparing Vincenty to alternatives
Karney's 2013 algorithm has displaced Vincenty in some modern geodesy software because it is guaranteed to converge for all input pairs and uses fewer trigonometric calls per pass. Vincenty remains the implementation in PostGIS, GeographicLib's legacy mode, and most aviation-flight-plan tools because it is well-understood, fast for typical inputs, and its 0.5 mm precision is two orders of magnitude finer than any GNSS receiver.
The latency table is a working-bench measurement, not a publication — it captures the ratio rather than the absolute numbers. Vincenty's per-call cost is roughly five times Haversine but lands within tens of microseconds of Karney for typical inputs. On the scale of a Next.js render cycle, none of these are the bottleneck.