Vincenty's formulae are two related iterative methods used in geodesy to calculate the distance between two points on the surface of a spheroid, developed by Thaddeus Vincenty (1975a). They are based on the assumption that the figure of the Earth is an oblate spheroid, and hence are more accurate than methods that assume a spherical Earth, such as great-circle distance. The first (direct) method computes the location of a point that is a given distance and azimuth (direction) from another point. The second (inverse) method computes the geographical distance and azimuth between two given points. They have been widely used in geodesy because they are accurate to within 0.5 mm (0.020 in) on the Earth ellipsoid.
Background Vincenty's goal was to express existing algorithms for geodesics on an ellipsoid in a form that minimized the program length (Vincenty 1975a). His unpublished report (1975b) mentions the use of a Wang 720 desk calculator, which had only a few kilobytes of memory. To obtain good accuracy for long lines, the solution uses the classical solution of Legendre (1806), Bessel (1825), and Helmert (1880) based on the auxiliary sphere. Vincenty relied on formulation of this method given by Rainsford, 1955. Legendre showed that an ellipsoidal geodesic can be exactly mapped to a great circle on the auxiliary sphere by mapping the geographic latitude to reduced latitude and setting the azimuth of the great circle equal to that of the geodesic. The longitude on the ellipsoid and the distance along the geodesic are then given in terms of the longitude on the sphere and the arc length along the great circle by simple integrals. Bessel and Helmert gave rapidly converging series for these integrals, which allow the geodesic to be computed with arbitrary accuracy. In order to minimize the program size, Vincenty took these series, re-expanded them using the first term of each series as the small parameter, and truncated them to O ( f 3 ) {\displaystyle O(f^{3})} . This resulted in compact expressions for the longitude and distance integrals. The expressions were put in Horner (or nested) form, since this allows polynomials to be evaluated using only a single temporary register. Finally, simple iterative techniques were used to solve the implicit equations in the direct and inverse methods; even though these are slow (and in the case of the inverse method it sometimes does not converge), they result in the least increase in code size.
Notation Define the following notation:
Inverse problem Given the coordinates of the two points (Φ1, L1) and (Φ2, L2), the inverse problem finds the azimuths α1, α2 and the ellipsoidal distance s. Calculate U1, U2 and L, and set initial value of λ = L. Then iteratively evaluate the following equations until λ converges:
sin σ = ( cos U 2 sin λ ) 2 + ( cos U 1 sin U 2 − sin U 1 cos U 2 cos λ ) 2 {\displaystyle \sin \sigma ={\sqrt {\left(\cos U_{2}\sin \lambda \right)^{2}+\left(\cos U_{1}\sin U_{2}-\sin U_{1}\cos U_{2}\cos \lambda \right)^{2}}}}
cos σ = sin U 1 sin U 2 + cos U 1 cos U 2 cos λ {\displaystyle \cos \sigma =\sin U_{1}\sin U_{2}+\cos U_{1}\cos U_{2}\cos \lambda \,}
σ = arctan2 ( sin σ , cos σ ) {\displaystyle \sigma =\operatorname {arctan2} \left(\sin \sigma ,\cos \sigma \right)}
sin α = cos U 1 cos U 2 sin λ sin σ {\displaystyle \sin \alpha ={\frac {\cos U_{1}\cos U_{2}\sin \lambda }{\sin \sigma }}}
cos 2 α = 1 − sin 2 α {\displaystyle \cos ^{2}\alpha =1-\sin ^{2}\alpha }
cos ( 2 σ m ) = cos σ − 2 sin U 1 sin U 2 cos 2 α = cos σ − 2 sin U 1 sin U 2 1 − sin 2 α {\displaystyle \cos \left(2\sigma _{\text{m}}\right)=\cos \sigma -{\frac {2\sin U_{1}\sin U_{2}}{\cos ^{2}\alpha }}=\cos \sigma -{\frac {2\sin U_{1}\sin U_{2}}{1-\sin ^{2}\alpha }}}
C = f 16 cos 2 α [ 4 + f ( 4 − 3 cos 2 α ) ] {\displaystyle C={\frac {f}{16}}\cos ^{2}\alpha \left[4+f\left(4-3\cos ^{2}\alpha \right)\right]}
λ = L + ( 1 − C ) f sin α { σ + C sin σ [ cos ( 2 σ m ) + C cos σ ( − 1 + 2 cos 2 ( 2 σ m ) ) ] } {\displaystyle \lambda =L+(1-C)f\sin \alpha \left\{\sigma +C\sin \sigma \left[\cos \left(2\sigma _{\text{m}}\right)+C\cos \sigma \left(-1+2\cos ^{2}\left(2\sigma _{\text{m}}\right)\right)\right]\right\}}
When λ has converged to the desired degree of accuracy (10−12 corresponds to approximately 0.006 mm), evaluate the following:
u 2 = cos 2 α ( a 2 − b 2 b 2 ) A = 1 + u 2 16384 ( 4096 + u 2 [ − 768 + u 2 ( 320 − 175 u 2 ) ] ) B = u 2 1024 ( 256 + u 2 [ − 128 + u 2 ( 74 − 47 u 2 ) ] ) Δ σ = B sin σ { cos ( 2 σ m ) + 1 4 B ( cos σ [ − 1 + 2 cos 2 ( 2 σ m ) ] − 1 6 B cos [ 2 σ m ] [ − 3 + 4 sin 2 σ ] [ − 3 + 4 cos 2 ( 2 σ m ) ] ) } s = b A ( σ − Δ σ ) α 1 = arctan2 ( cos U 2 sin λ , cos U 1 sin U 2 − sin U 1 cos U 2 cos λ ) α 2 = arctan2 ( cos U 1 sin λ , − sin U 1 cos U 2 + cos U 1 sin U 2 cos λ ) {\displaystyle {\begin{aligned}u^{2}&=\cos ^{2}\alpha \left({\frac {a^{2}-b^{2}}{b^{2}}}\right)\\A&=1+{\frac {u^{2}}{16384}}\left(4096+u^{2}\left[-768+u^{2}\left(320-175u^{2}\right)\right]\right)\\B&={\frac {u^{2}}{1024}}\left(256+u^{2}\left[-128+u^{2}\left(74-47u^{2}\right)\right]\right)\\\Delta \sigma &=B\sin \sigma \left\{\cos(2\sigma _{\text{m}})+{\frac {1}{4}}B\left(\cos \sigma \left[-1+2\cos ^{2}\left(2\sigma _{\text{m}}\right)\right]-{\frac {1}{6}}B\cos \left[2\sigma _{\text{m}}\right]\left[-3+4\sin ^{2}\sigma \right]\left[-3+4\cos ^{2}\left(2\sigma _{\text{m}}\right)\right]\right)\right\}\\s&=bA(\sigma -\Delta \sigma )\,\\\alpha _{1}&=\operatorname {arctan2} \left(\cos U_{2}\sin \lambda ,\cos U_{1}\sin U_{2}-\sin U_{1}\cos U_{2}\cos \lambda \right)\\\alpha _{2}&=\operatorname {arctan2} \left(\cos U_{1}\sin \lambda ,-\sin U_{1}\cos U_{2}+\cos U_{1}\sin U_{2}\cos \lambda \right)\end{aligned}}}
Between two nearly antipodal points, the iterative formula may fail to converge; this will occur when the first guess at λ as computed by the equation above is greater than π in absolute value.
Direct problem Given an initial point (Φ1, L1) and initial azimuth, α1, and a distance, s, along the geodesic the problem is to find the end point (Φ2, L2) and azimuth, α2. Start by calculating the following:
U 1 = arctan [ ( 1 − f ) tan ϕ 1 ] σ 1 = arctan2 ( tan U 1 , cos α 1 ) sin α = cos U 1 sin α 1 u 2 = cos 2 α ( a 2 − b 2 b 2 ) = ( 1 − sin 2 α ) ( a 2 − b 2 b 2 ) A = 1 + u 2 16384 ( 4096 + u 2 [ − 768 + u 2 ( 320 − 175 u 2 ) ] ) B = u 2 1024 ( 256 + u 2 [ − 128 + u 2 ( 74 − 47 u 2 ) ] ) {\displaystyle {\begin{aligned}U_{1}&=\arctan \left[(1-f)\tan \phi _{1}\right]\\\sigma _{1}&=\operatorname {arctan2} \left(\tan U_{1},\cos \alpha _{1}\right)\\\sin \alpha &=\cos U_{1}\sin \alpha _{1}\\u^{2}&=\cos ^{2}\alpha \left({\frac {a^{2}-b^{2}}{b^{2}}}\right)=\left(1-\sin ^{2}\alpha \right)\left({\frac {a^{2}-b^{2}}{b^{2}}}\right)\\A&=1+{\frac {u^{2}}{16384}}\left(4096+u^{2}\left[-768+u^{2}(320-175u^{2})\right]\right)\\B&={\frac {u^{2}}{1024}}\left(256+u^{2}\left[-128+u^{2}\left(74-47u^{2}\right)\right]\right)\end{aligned}}}
Then, using an initial value σ = s b A {\displaystyle \sigma ={\tfrac {s}{bA}}} , iterate the following equations until there is no significant change in σ:
2 σ m = 2 σ 1 + σ Δ σ = B sin σ { cos ( 2 σ m ) + 1 4 B ( cos σ [ − 1 + 2 cos 2 ( 2 σ m ) ] − 1 6 B cos [ 2 σ m ] [ − 3 + 4 sin 2 σ ] [ − 3 + 4 cos 2 ( 2 σ m ) ] ) } σ = s b A + Δ σ {\displaystyle {\begin{aligned}2\sigma _{\text{m}}&=2\sigma _{1}+\sigma \\\Delta \sigma &=B\sin \sigma \left\{\cos \left(2\sigma _{\text{m}}\right)+{\frac {1}{4}}B\left(\cos \sigma \left[-1+2\cos ^{2}\left(2\sigma _{\text{m}}\right)\right]-{\frac {1}{6}}B\cos \left[2\sigma _{\text{m}}\right]\left[-3+4\sin ^{2}\sigma \right]\left[-3+4\cos ^{2}\left(2\sigma _{\text{m}}\right)\right]\right)\right\}\\\sigma &={\frac {s}{bA}}+\Delta \sigma \end{aligned}}}
Once σ is obtained to sufficient accuracy evaluate:
ϕ 2 = arctan2 ( sin U 1 cos σ + cos U 1 sin σ cos α 1 , ( 1 − f ) sin 2 α + ( sin U 1 sin σ − cos U 1 cos σ cos α 1 ) 2 ) λ = arctan2 ( sin σ sin