The Vincenty Formula Explained

The great-circle-distance pillar introduces Vincenty's formula as the ellipsoidal upgrade to haversine. This article goes deeper: the math, the convergence behaviour, the antipodal edge case, and the Karney 2013 algorithm that's replacing Vincenty in modern libraries.

The need for an ellipsoidal formula

The spherical haversine formula assumes Earth is a sphere of some mean radius R. The real Earth (modelled as the WGS 84 ellipsoid) has a 21.4 km flattening at the poles, producing ~0.5% systematic error in spherical distance calculations.

For most use cases this is invisible. But:

• A 5,000 km transcontinental flight: ~25 km of error.
• A 10,000 km transoceanic flight: ~50 km of error.
• A 100 m surveying baseline: 0.5 m of error.

For surveying, aviation route optimisation, or any work that quotes distance to better than a percent, the spherical approximation is wrong. Vincenty's ellipsoidal formula removes this error.

The inverse problem

Vincenty's formula has two variants:

• Inverse: given two points (φ₁, λ₁) and (φ₂, λ₂), find the geodesic distance s between them and the initial and final azimuths.
• Direct: given a starting point (φ₁, λ₁), an initial azimuth α₁, and a distance s, find the destination point (φ₂, λ₂) and the final azimuth.

The inverse is what coordinate-distance computation uses; the direct is what dead-reckoning navigation uses. Vincenty published closed-form-with-iteration solutions for both.

The inverse iteration

The inverse formula starts by computing reduced latitudes U₁, U₂ (the latitudes on a sphere of the ellipsoid's polar radius), then iterates the difference in longitude on the auxiliary sphere λ until it converges. The full system of equations involves about a dozen intermediate variables computed each iteration; see Vincenty's original 1975 paper or any geodetic textbook for the complete derivation.

The iterative core, simplified:

1. Initialise λ = λ₂ − λ₁ (the longitude difference on the
   ellipsoid).
2. Compute sin(σ), cos(σ), σ (angular separation on auxiliary
   sphere); sin(α) (azimuth on auxiliary sphere); cos²(α);
   cos(2σₘ) (angular position of midpoint).
3. Compute C, the auxiliary integration constant.
4. Update λ using the iteration formula:
   λ_new = (λ₂ − λ₁) + (1 − C)·f·sin(α) · (σ + C·sin(σ) · (cos(2σₘ) + C·cos(σ) · (−1 + 2·cos²(2σₘ))))
5. If |λ_new − λ| > tolerance (typically 10⁻¹²), set λ = λ_new
   and go to step 2.
6. When converged, compute the auxiliary u², A, B, and finally
   the geodesic distance:
   s = b · A · (σ − Δσ)
   where b is the polar radius and Δσ is the high-order correction.

The iteration converges to roughly 10⁻¹² radians of angular precision in 5–10 cycles for most input pairs. Each iteration involves about 5–10 trig calls plus standard arithmetic; total cost is around a few microseconds in a typical implementation.

The Coordinately library src/lib/coords/vincenty.ts is a TypeScript port of Vincenty's original Algol-60 code with the standard convergence-check + haversine-fallback logic.

A worked example

The Empire State Building (40.7484°N, 73.9857°W) to the Statue of Liberty (40.6892°N, 74.0445°W). Running Vincenty on WGS 84:

Inputs:
  φ₁ = 40.7484°,  λ₁ = −73.9857°
  φ₂ = 40.6892°,  λ₂ = −74.0445°
  a = 6,378,137.0 m,  f = 1/298.257223563

Iteration converges in ~5 cycles to:
  σ = 7.6 × 10⁻⁴ rad
  α = 230.6° (initial azimuth from Empire State toward Statue of Liberty)
  α' = 230.6° (final azimuth on arrival)

Geodesic distance:
  s = 4,838 m

Compare to haversine: 4,832 m. The 6 m discrepancy is the sphere-vs-ellipsoid difference for this 5 km southwest path.

For the longer JFK → LHR case, the discrepancy grows to ~46 m (Vincenty 5,539 km vs haversine 5,572 km — note the haversine overestimates here because the path crosses high latitudes where the spherical assumption overstates path length).

The antipodal failure mode

Vincenty's well-known limitation: when the two endpoints are within about 50 km of being exactly antipodal (geometrically opposite on Earth, ~20,000 km apart), the iteration fails to converge. The geometric reason: at antipodes, every great-circle- like path between the points is the same length; the iteration has no unique direction to converge toward, and the longitude update oscillates instead of settling.

The standard remedy:

1. Track the change in λ between iterations.
2. If the iteration exceeds some maximum (typically 100–200 iterations) without converging, declare non-convergence.
3. Fall back to the spherical haversine formula, which always produces a result.

The fall-back result has ~0.5% accuracy rather than millimetre accuracy — acceptable because anti-podal-pair distance is rarely needed at sub-percent precision anyway (no one is engineering a project that crosses Earth's antipode).

The Coordinately tools detect Vincenty non-convergence after 200 iterations and emit the haversine result with a method-name flag in the response metadata.

Karney 2013 — the modern successor

Charles Karney's 2013 paper Algorithms for geodesics publishes a robust ellipsoidal-geodesic algorithm that converges everywhere — including at antipodes. The algorithm is algorithmically more sophisticated than Vincenty: it uses series expansions in different variables and handles the antipodal case through reformulation rather than fallback.

Properties:

• Accuracy: nanometre precision (much better than Vincenty's millimetre).
• Robustness: converges for any pair of points on the ellipsoid; no antipodal failure mode.
• Performance: comparable to Vincenty in iteration count; similar per-call cost in typical implementations.
• Implementation: more complex than Vincenty; the GeographicLib library is the reference implementation.

Modern GIS (PROJ 6+, PostGIS via PROJ, GDAL, QGIS) increasingly use Karney through GeographicLib. The Coordinately codebase currently uses Vincenty with haversine fallback because it's the simpler implementation and the antipodal case is rare in practice; a future migration to Karney is on the roadmap.

A short note about Vincenty the person. Thaddeus Vincenty (1920–2002) was a Polish-American geodesist who worked at the US Defense Mapping Agency (later NGA) and contributed several foundational papers on coordinate computation, of which the 1975 geodesic formula is the most-cited. The formula was originally implemented in Algol 60 in his paper's appendix; modern implementations are direct translations to C, Fortran, Python, JavaScript, and other languages. The mathematics has remained unchanged since publication; the implementations have proliferated into every major geospatial library.

The direct problem

Vincenty's direct formula is the inverse of the inverse: given a starting point, an azimuth, and a distance, find the endpoint. The iteration structure is similar but the variables are different — the auxiliary sphere's angular separation σ is iterated instead of the longitude difference λ.

The direct formula is what dead-reckoning navigation uses: given a known starting point and the navigator's recorded course and distance, project the current position. It's also used in waypoint-projection routines (find the point 100 km west- northwest of New York, where “west-northwest” is an azimuth).

The Coordinately src/lib/coords/vincenty.ts implements both inverse and direct.

Comparison with other formulas

| Formula | Earth model | Accuracy | Robust everywhere? | Speed | | ------------------------- | ---------------- | ------------------- | ------------------ | ------------- | | Spherical law of cosines | Sphere | ~0.5% (unstable short) | Yes | Fast (single call) | | Haversine | Sphere | ~0.5% (stable) | Yes | Fast (single call) | | Vincenty inverse | WGS 84 | ~1 mm | Fails near antipodes | Iterative (5–10) | | Karney 2013 | WGS 84 | ~1 nm | Yes | Iterative (similar) |

For sub-metre work where antipodal robustness matters, Karney is the best choice. For sub-metre work where antipodes are guaranteed not to occur, Vincenty with no fallback is fine. For ~0.5% accuracy at low computational cost, haversine is the right choice.

Implementation pitfalls

A few common Vincenty-implementation mistakes:

Forgetting the convergence cap. An implementation without an iteration limit hangs forever on antipodal inputs. Always cap the iteration at ~100–200 cycles and treat non-convergence as a signal to fall back.

Wrong ellipsoid parameters. Vincenty needs (a, f). Using WGS 84 (a = 6,378,137.0 m, f = 1/298.257223563) for one ellipsoid and GRS80 for another produces sub-millimetre discrepancies — usually negligible, but worth being aware of.

Returning negative distances or NaN. The iteration can produce edge cases (input points equal, points exactly on a pole, points across the antimeridian) where defensive checks are needed. Coordinately's implementation has explicit handling for each.

Mixing degrees and radians. As with haversine, all trig calls must be in radians; latitudes / longitudes typically come in degrees. Convert at the boundary.

Recovering azimuth from atan2(y, x). The azimuth output is in [−π, π]; normalise to [0, 2π] or to degree-compass convention as needed before display.

Common misconceptions

“Vincenty is the most accurate.” Karney 2013 is more accurate (nanometres vs millimetres). For practical work the difference doesn't matter, but for academic citation accuracy claims, Karney is the modern reference.

“Vincenty's antipodal failure is a deal-breaker.” Only if you compute distances between antipodes. In practice, the antipodal failure affects an extremely narrow band of input pairs that real applications rarely encounter.

“Vincenty needs the geoid.” It needs the ellipsoid (via a and f) — not the geoid. Vincenty operates on the geometric ellipsoid surface, not the equipotential geoid. Vertical considerations (orthometric vs ellipsoidal height) are separate.

“Vincenty distances change with elevation.” The formula computes distance along the ellipsoid surface, ignoring elevation. For a flight at 35,000 feet, the actual travelled distance is slightly longer than the Vincenty surface distance (by roughly the geocentric-radius effect, ~0.2% at typical cruise altitude). For ground-based surveying this is negligible.

“Vincenty is slow.” A few microseconds per call. Roughly 10× slower than haversine in raw compute terms, but in the context of a typical web request (database lookup, network I/O, rendering) the cost is invisible.

“Vincenty is deprecated.” Not quite. Karney is the modern academic reference, but Vincenty remains the most- implemented ellipsoidal geodesic formula in the world. It's in every GIS, every navigation software, every coordinate library. Modern usage trends toward Karney for new implementations; Vincenty is the legacy default and the standard reference.

“Vincenty needs special hardware.” No — it's pure floating-point arithmetic that runs on any CPU with a math library. Vincenty's original 1975 implementation ran on a mainframe in Algol 60 with no special facilities; modern implementations are sub-microsecond per call on commodity hardware. The historical record of Vincenty's formula producing centimetre-accurate results in pre-GPS-era surveying is a testament to its mathematical robustness, not its hardware requirements.