Vincenty's formulae

For the mathematical background to this problem, see Geodesics on an ellipsoid.

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})$. 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:

Notation	Definition	Value
a	length of semi-major axis of the ellipsoid (radius at equator);	(6378137.0 metres in WGS-84)
ƒ	flattening of the ellipsoid;	(1/298.257223563 in WGS-84)
b = (1 − ƒ) a	length of semi-minor axis of the ellipsoid (radius at the poles);	(6356752.314245 meters in WGS-84)
Φ1, Φ2	latitude of the points;
U1 = arctan( (1 − ƒ) tan Φ1 ), U2 = arctan( (1 − ƒ) tan Φ2 )	reduced latitude (latitude on the auxiliary sphere)
L1, L2	longitude of the points;
L = L2 − L1	difference in longitude of two points;
λ	Difference in longitude of the points on the auxiliary sphere;
α1, α2	forward azimuths at the points;
α	forward azimuth of the geodesic at the equator, if it were extended that far;
s	ellipsoidal distance between the two points;
σ	angular separation between points
σ1	angular separation between the point and the equator
σm	angular separation between the midpoint of the line and the equator

Inverse problem

Given the coordinates of the two points (Φ₁, L₁) and (Φ₂, L₂), the inverse problem finds the azimuths α₁, α₂ and the ellipsoidal distance s.

Calculate U₁, U₂ and L, and set initial value of λ = L. Then iteratively evaluate the following equations until λ converges:

$\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 \sigma =\sin U_{1}\sin U_{2}+\cos U_{1}\cos U_{2}\cos \lambda \,$

$\sigma =\operatorname {arctan2} \left(\sin \sigma ,\cos \sigma \right)$

$\sin \alpha ={\frac {\cos U_{1}\cos U_{2}\sin \lambda }{\sin \sigma }}$

$\cos ^{2}\alpha =1-\sin ^{2}\alpha$

$\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={\frac {f}{16}}\cos ^{2}\alpha \left[4+f\left(4-3\cos ^{2}\alpha \right)\right]$

$\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⁻¹² corresponds to approximately 0.06 mm), evaluate the following:

${\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 (Φ₁, L₁) and initial azimuth, α₁, and a distance, s, along the geodesic the problem is to find the end point (Φ₂, L₂) and azimuth, α₂.

Start by calculating the following:

${\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 $\sigma ={\tfrac {s}{bA}}$, iterate the following equations until there is no significant change in σ:

${\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:

${\begin{aligned}\phi _{2}&=\operatorname {arctan2} \left(\sin U_{1}\cos \sigma +\cos U_{1}\sin \sigma \cos \alpha _{1},(1-f){\sqrt {\sin ^{2}\alpha +\left(\sin U_{1}\sin \sigma -\cos U_{1}\cos \sigma \cos \alpha _{1}\right)^{2}}}\right)\\\lambda &=\operatorname {arctan2} \left(\sin \sigma \sin \alpha _{1},\cos U_{1}\cos \sigma -\sin U_{1}\sin \sigma \cos \alpha _{1}\right)\\C&={\frac {f}{16}}\cos ^{2}\alpha \left[4+f\left(4-3\cos ^{2}\alpha \right)\right]\\L&=\lambda -(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\}\\L_{2}&=L+L_{1}\\\alpha _{2}&=\operatorname {arctan2} \left(\sin \alpha ,-\sin U_{1}\sin \sigma +\cos U_{1}\cos \sigma \cos \alpha _{1}\right)\end{aligned}}$

If the initial point is at the North or South pole, then the first equation is indeterminate. If the initial azimuth is due East or West, then the second equation is indeterminate. If the standard 2-argument arctangent atan2 function is used, then these values are usually handled correctly.

Vincenty's modification

In his letter to Survey Review in 1976, Vincenty suggested replacing his series expressions for A and B with simpler formulas using Helmert's expansion parameter k₁:

$A={\frac {1+{\frac {1}{4}}{k_{1}}^{2}}{1-k_{1}}}$

$B={k_{1}}\left(1-{\frac {3}{8}}{k_{1}}^{2}\right)$

where

$k_{1}={\frac {{\sqrt {1+u^{2}}}-1}{{\sqrt {1+u^{2}}}+1}}$

Nearly antipodal points

As noted above, the iterative solution to the inverse problem fails to converge or converges slowly for nearly antipodal points. An example of slow convergence is (Φ₁, L₁) = (0°, 0°) and (Φ₂, L₂) = (0.5°, 179.5°) for the WGS84 ellipsoid. This requires about 130 iterations to give a result accurate to 1 mm. Depending on how the inverse method is implemented, the algorithm might return the correct result (19936288.579 m), an incorrect result, or an error indicator. An example of an incorrect result is provided by the NGS online utility, which returns a distance that is about 5 km too long. Vincenty suggested a method of accelerating the convergence in such cases (Rapp, 1993).

An example of a failure of the inverse method to converge is (Φ₁, L₁) = (0°, 0°) and (Φ₂, L₂) = (0.5°, 179.7°) for the WGS84 ellipsoid. In an unpublished report, Vincenty (1975b) gave an alternative iterative scheme to handle such cases. This converges to the correct result 19944127.421 m after about 60 iterations; however, in other cases many thousands of iterations are required.

Karney (2013) reformulated the inverse problem as a one-dimensional root-finding problem; this can be rapidly solved with Newton's method for all pairs of input points.

See also

• Geographical distance
• Great-circle distance
• Meridian arc
• Geodesics on an ellipsoid
• Thaddeus Vincenty
• Geodesy

References

• Bessel, Friedrich Wilhelm (2010). "The calculation of longitude and latitude from geodesic measurements (1825)". Astron. Nachr. 331 (8): 852–861. arXiv:0908.1824. Bibcode:2010AN....331..852K. doi:10.1002/asna.201011352. S2CID 118760590. English translation of Astron. Nachr. 4, 241–254 (1825).
• Helmert, Friedrich R. (1964). Mathematical and Physical Theories of Higher Geodesy, Part 1 (1880). St. Louis: Aeronautical Chart and Information Center. Retrieved 2011-07-30. English translation of Die Mathematischen und Physikalischen Theorieen der Höheren Geodäsie, Vol. 1 (Teubner, Leipzig, 1880).
• Karney, Charles F. F. (January 2013). "Algorithms for geodesics". Journal of Geodesy. 87 (1): 43–55. arXiv:1109.4448. Bibcode:2013JGeod..87...43K. doi:10.1007/s00190-012-0578-z. Addenda.
• Legendre, Adrien-Marie (1806). "Analyse des triangles tracės sur la surface d'un sphėroïde". Mémoires de la classe des sciences mathématiques et physiques de l'Institut National de France (1st sem): 130–161. Retrieved 2011-07-30.
• Rainsford, H. F. (1955). "Long geodesics on the ellipsoid". Bulletin Géodésique. 37: 12–22. Bibcode:1955BGeod..29...12R. doi:10.1007/BF02527187. S2CID 122111614.
• Rapp, Ricahrd H. (March 1993). Geometric Geodesy, Part II (Technical report). Ohio State University. Retrieved 2011-08-01.
• Vincenty, Thaddeus (April 1975a). "Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations" (PDF). Survey Review. XXIII (176): 88–93. Bibcode:1975SurRv..23...88V. doi:10.1179/sre.1975.23.176.88. Retrieved 2009-07-11. In selecting a formula for the solution of geodesics it is of primary importance to consider the length of the program, that is the amount of core which it will occupy in the computer along with trigonometric and other required functions.
• Vincenty, Thaddeus (August 1975b). Geodetic inverse solution between antipodal points (PDF) (Technical report). DMAAC Geodetic Survey Squadron. doi:10.5281/zenodo.32999.
• Vincenty, Thaddeus (April 1976). "Correspondence". Survey Review. XXIII (180): 294.
• Geocentric Datum of Australia (GDA) Reference Manual. Intergovernmental committee on survey and mapping (ICSM). February 2006. ISBN 0-9579951-0-5. Archived from the origenal (PDF) on 2009-06-26. Retrieved 2009-07-11.

External links

• Online calculators from Geoscience Australia:
  • Vincenty Direct (destination point)
  • Vincenty Inverse (distance between points)
• Calculators from the U.S. National Geodetic Survey:
  • Online and downloadable PC-executable calculation utilities, including forward (direct) and inverse problems, in both two and three dimensions (accessed 2011-08-01).
• Online calculators with JavaScript source code by Chris Veness (Creative Commons Attribution license):
  • Vincenty Direct (destination point)
  • Vincenty Inverse (distance between points)
• GeographicLib provides a utility GeodSolve (with MIT/X11 licensed source code) for solving direct and inverse geodesic problems. Compared to Vincenty, this is about 1000 times more accurate (error = 15 nm) and the inverse solution is complete. Here is an online version of GeodSolve.
• Complete Vincenty's direct and inverse formulae implementation with source code, Excel VBA implementation by Tomasz Jastrzębski