Vincenty's formulae

Last updated September 21, 2024

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,^{[ clarification needed ]} 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)
Φ₁, Φ₂	latitude of the points;
$U 1 = arctan( (1 - ƒ) tan Φ 1), U 2 = arctan( (1 - ƒ) tan Φ 2)$	reduced latitude (latitude on the auxiliary sphere)
L₁, L₂	longitude of the points;
$L = L 2 - L 1$	difference in longitude of two points;
λ	Difference in longitude of the points on the auxiliary sphere;
α₁, α₂	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
σ₁	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)

^[1]

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

^[2]

\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 }}

^[3]

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.^{[ clarification needed ]}

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.

Notes

↑ σ is not evaluated directly from sin σ or cos σ to preserve numerical accuracy near the poles and equator
↑ If sin σ = 0 the value of sin α is indeterminate. It represents an end point coincident with, or diametrically opposed to, the start point.
↑ Where the start and end point are on the equator, C = 0 and the value of $\cos \left(2\sigma _{\text{m}}\right)$ is not used. The limiting value is $\cos \left(2\sigma _{\text{m}}\right)=-1$ .

Related Research Articles

In mathematical physics and mathematics, the Pauli matrices are a set of three $2 \times 2$ complex matrices that are traceless, Hermitian, involutory and unitary. Usually indicated by the Greek letter sigma, they are occasionally denoted by tau when used in connection with isospin symmetries.

An azimuth is the horizontal angle from a cardinal direction, most commonly north, in a local or observer-centric spherical coordinate system.

In physics and mathematics, the Lorentz group is the group of all Lorentz transformations of Minkowski spacetime, the classical and quantum setting for all (non-gravitational) physical phenomena. The Lorentz group is named for the Dutch physicist Hendrik Lorentz.

In mathematics, a Gaussian function, often simply referred to as a Gaussian, is a function of the base form $and with parametric extension for arbitrary real constants a, b and non-zero c . It is named after the mathematician Carl Friedrich Gauss. The graph of a Gaussian is a characteristic symmetric "bell curve" shape. The parameter a is the height of the curve's peak, b is the position of the center of the peak, and c controls the width of the "bell".$

The great-circle distance, orthodromic distance, or spherical distance is the distance between two points on a sphere, measured along the great-circle arc between them. This arc is the shortest path between the two points on the surface of the sphere.

In geodesy, conversion among different geographic coordinate systems is made necessary by the different geographic coordinate systems in use across the world and over time. Coordinate conversion is composed of a number of different types of conversion: format change of geographic coordinates, conversion of coordinate systems, or transformation to different geodetic datums. Geographic coordinate conversion has applications in cartography, surveying, navigation and geographic information systems.

In particle physics, the Klein–Nishina formula gives the differential cross section of photons scattered from a single free electron, calculated in the lowest order of quantum electrodynamics. It was first derived in 1928 by Oskar Klein and Yoshio Nishina, constituting one of the first successful applications of the Dirac equation. The formula describes both the Thomson scattering of low energy photons and the Compton scattering of high energy photons, showing that the total cross section and expected deflection angle decrease with increasing photon energy.

<span class="mw-page-title-main">Universal Transverse Mercator coordinate system</span> Map projection system

The Universal Transverse Mercator (UTM) is a map projection system for assigning coordinates to locations on the surface of the Earth. Like the traditional method of latitude and longitude, it is a horizontal position representation, which means it ignores altitude and treats the earth surface as a perfect ellipsoid. However, it differs from global latitude/longitude in that it divides earth into 60 zones and projects each to the plane as a basis for its coordinates. Specifying a location means specifying the zone and the x, y coordinate in that plane. The projection from spheroid to a UTM zone is some parameterization of the transverse Mercator projection. The parameters vary by nation or region or mapping system.

<span class="mw-page-title-main">Great-circle navigation</span> Flight or sailing route along the shortest path between two points on a globes surface

Great-circle navigation or orthodromic navigation is the practice of navigating a vessel along a great circle. Such routes yield the shortest distance between two points on the globe.

Ellipsoidal coordinates are a three-dimensional orthogonal coordinate system $that generalizes the two-dimensional elliptic coordinate system. Unlike most three-dimensional orthogonal coordinate systems that feature quadratic coordinate surfaces, the ellipsoidal coordinate system is based on confocal quadrics.$

In probability theory, the inverse Gaussian distribution is a two-parameter family of continuous probability distributions with support on (0,∞).

The Jaynes–Cummings model is a theoretical model in quantum optics. It describes the system of a two-level atom interacting with a quantized mode of an optical cavity, with or without the presence of light. It was originally developed to study the interaction of atoms with the quantized electromagnetic field in order to investigate the phenomena of spontaneous emission and absorption of photons in a cavity.

In computational chemistry, a constraint algorithm is a method for satisfying the Newtonian motion of a rigid body which consists of mass points. A restraint algorithm is used to ensure that the distance between mass points is maintained. The general steps involved are: (i) choose novel unconstrained coordinates, (ii) introduce explicit constraint forces, (iii) minimize constraint forces implicitly by the technique of Lagrange multipliers or projection methods.

A ratio distribution is a probability distribution constructed as the distribution of the ratio of random variables having two other known distributions. Given two random variables X and Y, the distribution of the random variable Z that is formed as the ratio Z = X/Y is a ratio distribution.

<span class="mw-page-title-main">Great ellipse</span> Ellipse on a spheroid centered on its origin

A great ellipse is an ellipse passing through two points on a spheroid and having the same center as that of the spheroid. Equivalently, it is an ellipse on the surface of a spheroid and centered on the origin, or the curve formed by intersecting the spheroid by a plane through its center. For points that are separated by less than about a quarter of the circumference of the earth, about $, the length of the great ellipse connecting the points is close to the geodesic distance. The great ellipse therefore is sometimes proposed as a suitable route for marine navigation. The great ellipse is special case of an earth section path.$

In mathematics, the spectral theory of ordinary differential equations is the part of spectral theory concerned with the determination of the spectrum and eigenfunction expansion associated with a linear ordinary differential equation. In his dissertation, Hermann Weyl generalized the classical Sturm–Liouville theory on a finite closed interval to second order differential operators with singularities at the endpoints of the interval, possibly semi-infinite or infinite. Unlike the classical case, the spectrum may no longer consist of just a countable set of eigenvalues, but may also contain a continuous part. In this case the eigenfunction expansion involves an integral over the continuous part with respect to a spectral measure, given by the Titchmarsh–Kodaira formula. The theory was put in its final simplified form for singular differential equations of even degree by Kodaira and others, using von Neumann's spectral theorem. It has had important applications in quantum mechanics, operator theory and harmonic analysis on semisimple Lie groups.

<span class="mw-page-title-main">Normal-inverse-gamma distribution</span>

In probability theory and statistics, the normal-inverse-gamma distribution is a four-parameter family of multivariate continuous probability distributions. It is the conjugate prior of a normal distribution with unknown mean and variance.

<span class="mw-page-title-main">Geographical distance</span> Distance measured along the surface of the Earth

Geographical distance or geodetic distance is the distance measured along the surface of the Earth, or the shortest arch length.

<span class="mw-page-title-main">Geodesics on an ellipsoid</span> Shortest paths on a bounded deformed sphere-like quadric surface

The study of geodesics on an ellipsoid arose in connection with geodesy specifically with the solution of triangulation networks. The figure of the Earth is well approximated by an oblate ellipsoid, a slightly flattened sphere. A geodesic is the shortest path between two points on a curved surface, analogous to a straight line on a plane surface. The solution of a triangulation network on an ellipsoid is therefore a set of exercises in spheroidal trigonometry.

<span class="mw-page-title-main">Earth section paths</span> Plane curved by the intersection of an earth ellipsoid and a plane

Earth section paths are plane curves defined by the intersection of an earth ellipsoid and a plane. Common examples include the great ellipse and normal sections. Earth section paths are useful as approximate solutions for geodetic problems, the direct and inverse calculation of geographic distances. The rigorous solution of geodetic problems involves skew curves known as geodesics.

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 original (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

This page is based on this Wikipedia article
Text is available under the CC BY-SA 4.0 license; additional terms may apply.
Images, videos and audio are available under their respective licenses.

[1] σ is not evaluated directly from sin σ or cos σ to preserve numerical accuracy near the poles and equator

[opposite-2] If sin σ = 0 the value of sin α is indeterminate. It represents an end point coincident with, or diametrically opposed to, the start point.

[equator-3] Where the start and end point are on the equator, C = 0 and the value of $\cos \left(2\sigma _{\text{m}}\right)$ is not used. The limiting value is $\cos \left(2\sigma _{\text{m}}\right)=-1$ .

[1]

[2]

[3]