Vincenty Inverse Method — Fully Worked Example
1) Input coordinates
Start: Les Sables-d’Olonne
φ₁ = 46.494953° → 0.81149001541 rad
λ₁ = −1.792091° → -0.03128378623 rad
End: Saint-François
φ₂ = 16.252360° → 0.28365719322 rad
λ₂ = −61.273320° → -1.06942707541 rad
$$L = \lambda_2 - \lambda_1 = -1.06942707541 - (-0.03128378623) = -1.03814328918 \text{ rad}$$
L = initial difference in longitude (radians)
2) WGS-84 ellipsoid parameters
$$a = 6378137\ \text{m}, \quad f = \frac{1}{298.257223563}, \quad b = a(1-f) = 6356752.314245\ \text{m}$$
a = equatorial radius, b = polar radius, f = flattening
3) Reduced latitudes
$$U_1 = \arctan((1-f) \tan\varphi_1) = \arctan(0.996647189335 \cdot \tan(0.81149001541)) = 0.80981293556 \text{ rad}$$
U₁ = latitude of start projected onto auxiliary sphere
$$U_2 = \arctan((1-f) \tan\varphi_2) = \arctan(0.996647189335 \cdot \tan(0.28365719322)) = 0.28275610843 \text{ rad}$$
U₂ = latitude of end projected onto auxiliary sphere
4) Iterative solution for λ
$$\lambda_0 = L = -1.03814328918$$
Initial guess for λ
$$\sin\sigma = \sqrt{(\cos U_2 \sin \lambda)^2 + (\cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos \lambda)^2}$$
σ = angular distance between points on auxiliary sphere
$$\cos\sigma = \sin U_1 \sin U_2 + \cos U_1 \cos U_2 \cos \lambda$$
cosσ = cosine of angular distance
$$\sigma = \arctan2(\sin\sigma, \cos\sigma)$$
σ = spherical arc length (radians)
$$\sin\alpha = \frac{\cos U_1 \cos U_2 \sin\lambda}{\sin\sigma}, \quad \cos^2\alpha = 1 - \sin^2\alpha$$
α = azimuth of geodesic; cos²α = squared north-south component
$$\cos2\sigma_m = \cos\sigma - \frac{2 \sin U_1 \sin U_2}{\cos^2\alpha}$$
σₘ = midpoint angular distance, used for ellipsoid corrections
$$C = \frac{f}{16} \cos^2\alpha (4 + f(4 - 3 \cos^2\alpha))$$
C = small coefficient for λ correction
$$\lambda_{i+1} = L + (1-C)f \sin\alpha \left[ \sigma + C \sin\sigma (\cos2\sigma_m + C \cos\sigma(-1 + 2 \cos^2 2\sigma_m)) \right]$$
λ is updated iteratively until Δλ < 10⁻¹² rad
| Iter | λ (rad) | Δλ | sinσ | cosσ | σ | sinα | cos²α | cos2σₘ |
|---|
| 0 | -1.03814328918 | — | — | — | — | — | — | — |
| 1 | -1.04041711352 | -0.00227382433 | 0.839510 | 0.543310 | 0.993219 | -0.672447 | 0.547387 | -0.217723 |
| 2 | -1.04042141420 | -0.00000430069 | 0.843102 | 0.537903 | 1.001753 | -0.676899 | 0.542038 | -0.210051 |
| 3 | -1.04042142233 | -8.13e-9 | 0.843553 | 0.537045 | 1.003865 | -0.677215 | 0.541379 | -0.209354 |
| 4 | -1.04042142235 | -1.54e-11 | 0.843553 | 0.537046 | 1.003866 | -0.677215 | 0.541379 | -0.209354 |
| 5 | -1.04042142235 | -2.91e-14 | 0.8435532581 | 0.5370455295 | 1.00386554952 | -0.677215389 | 0.541379317 | -0.209353772 |
Δλ = λᵢ - λᵢ₋₁; iteration stops when Δλ is very small
5) Ellipsoidal distance calculation
$$u^2 = \cos^2\alpha \frac{a^2 - b^2}{b^2} = 0.00364862414$$
u² = flattening factor effect
$$A = 1 + \frac{u^2}{16384}(4096 + u^2(-768 + u^2(320 - 175 u^2))) = 1.000911533$$
A = scale factor for distance
$$B = \frac{u^2}{1024}(256 + u^2(-128 + u^2(74 - 47 u^2))) = 0.0009104955$$
B = correction factor for Δσ
$$\Delta\sigma = B \sin\sigma \left[\cos2\sigma_m + \frac{B}{4} (\cos\sigma(-1 + 2\cos^2 2\sigma_m) - \frac{B}{6} \cos2\sigma_m(-3 + 4\sin^2\sigma)(-3 + 4\cos^2 2\sigma_m)) \right] = -0.00016088012$$
Δσ = additional angular distance correction due to flattening
$$s = b A (\sigma - \Delta\sigma) = 6356752.314245 \cdot 1.000911533 \cdot (1.00386554952 - (-0.00016088012)) = 6,388,165.05\ \text{m} = 6,388.17\ \text{km}$$
s = final geodesic distance along ellipsoid
$$s_{NM} = \frac{s}{1852} = 3,449.33\ \text{NM}$$
6) Initial and final bearings
$$\alpha_1 = \arctan2(\cos U_2 \sin\lambda, \cos U_1 \sin U_2 - \sin U_1 \cos U_2 \cos\lambda) = 259.11^\circ$$
α₁ = initial course at start point
$$\alpha_2 = \arctan2(\cos U_1 \sin\lambda, -\sin U_1 \cos U_2 + \cos U_1 \sin U_2 \cos\lambda) = 224.85^\circ$$
α₂ = final course at end point
7) Final summary
- Distance: 6,388 km
- Distance: 3,449 NM
- Initial bearing: 259.11°
- Final bearing: 224.85°
- High accuracy: millimeter level (non-antipodal)