Jump to

### Introduction

Consider a ellipsoid of revolution with equatorial radius *a*, polar
semi-axis *b*, and flattening *f* = (*a* − *b*)/*a* . Points on
the surface of the ellipsoid are characterized by their latitude φ
and longitude λ. (Note that latitude here means the
*geographical latitude*, the angle between the normal to the ellipsoid
and the equatorial plane).

The shortest path between two points on the ellipsoid at
(φ_{1}, λ_{1}) and (φ_{2},
λ_{2}) is called the geodesic. Its length is
*s*_{12} and the geodesic from point 1 to point 2 has forward
azimuths α_{1} and α_{2} at the two end
points. In this figure, we have λ_{12} =
λ_{2} − λ_{1}.

Traditionally two geodesic problems are considered:

the direct problem — given φ

_{1}, λ_{1}, α_{1},*s*_{12}, determine φ_{2}, λ_{2}, and α_{2}; this is solved by Geodesic.Direct.the inverse problem — given φ

_{1}, λ_{1}, φ_{2}, λ_{2}, determine*s*_{12}, α_{1}, and α_{2}; this is solved by Geodesic.Inverse.

### Additional properties

The routines also calculate several other quantities of interest

*S*_{12}is the area between the geodesic from point 1 to point 2 and the equator; i.e., it is the area, measured counter-clockwise, of the quadrilateral with corners (φ_{1},λ_{1}), (0,λ_{1}), (0,λ_{2}), and (φ_{2},λ_{2}). It is given in meters^{2}.*m*_{12}, the reduced length of the geodesic is defined such that if the initial azimuth is perturbed by*d*α_{1}(radians) then the second point is displaced by*m*_{12}*d*α_{1}in the direction perpendicular to the geodesic.*m*_{12}is given in meters. On a curved surface the reduced length obeys a symmetry relation,*m*_{12}+*m*_{21}= 0. On a flat surface, we have*m*_{12}=*s*_{12}.*M*_{12}and*M*_{21}are geodesic scales. If two geodesics are parallel at point 1 and separated by a small distance*dt*, then they are separated by a distance*M*_{12}*dt*at point 2.*M*_{21}is defined similarly (with the geodesics being parallel to one another at point 2).*M*_{12}and*M*_{21}are dimensionless quantities. On a flat surface, we have*M*_{12}=*M*_{21}= 1.- σ
_{12}is the arc length on the auxiliary sphere. This is a construct for converting the problem to one in spherical trigonometry. The spherical arc length from one equator crossing to the next is always 180°.

If points 1, 2, and 3 lie on a single geodesic, then the following addition rules hold:

*s*_{13}=*s*_{12}+*s*_{23}- σ
_{13}= σ_{12}+ σ_{23} *S*_{13}=*S*_{12}+*S*_{23}*m*_{13}=*m*_{12}*M*_{23}+*m*_{23}*M*_{21}*M*_{13}=*M*_{12}*M*_{23}− (1 −*M*_{12}*M*_{21})*m*_{23}/*m*_{12}*M*_{31}=*M*_{32}*M*_{21}− (1 −*M*_{23}*M*_{32})*m*_{12}/*m*_{23}

### Multiple shortest geodesics

The shortest distance found by solving the inverse problem is (obviously) uniquely defined. However, in a few special cases there are multiple azimuths which yield the same shortest distance. Here is a catalog of those cases:

- φ
_{1}= −φ_{2}(with neither point at a pole). If α_{1}= α_{2}, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting [α_{1},α_{2}] ← [α_{2},α_{1}], [*M*_{12},*M*_{21}] ← [*M*_{21},*M*_{12}],*S*_{12}← −*S*_{12}. (This occurs when the longitude difference is near ±180° for oblate ellipsoids.) - λ
_{2}= λ_{1}± 180° (with neither point at a pole). If α_{1}= 0° or ±180°, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting [α_{1},α_{2}] ← [−α_{1},−α_{2}],*S*_{12}← −*S*_{12}. (This occurs when φ_{2}is near −φ_{1}for prolate ellipsoids.) - Points 1 and 2 at opposite poles. There are infinitely many
geodesics which can be generated by setting
[α
_{1},α_{2}] ← [α_{1},α_{2}] + [δ,−δ], for arbitrary δ. (For spheres, this prescription applies when points 1 and 2 are antipodal.) *s*_{12}= 0 (coincident points). There are infinitely many geodesics which can be generated by setting [α_{1},α_{2}] ← [α_{1},α_{2}] + [δ,δ], for arbitrary δ.

### Background

The algorithms implemented by this package are given in Karney (2013) and are based on Bessel (1825) and Helmert (1880); the algorithm for areas is based on Danielsen (1989). These improve on the work of Vincenty (1975) in the following respects:

- The results are accurate to round-off for terrestrial ellipsoids (the error in the distance is less then 15 nanometers, compared to 0.1 mm for Vincenty).
- The solution of the inverse problem is always found. (Vincenty's method fails to converge for nearly antipodal points.)
- The routines calculate differential and integral properties of a geodesic. This allows, for example, the area of a geodesic polygon to be computed.

### References

- F. W. Bessel,
The calculation of longitude and latitude from geodesic measurements (1825),
Astron. Nachr.
**331**(8), 852–861 (2010), translated by C. F. F. Karney and R. E. Deakin. - F. R. Helmert, Mathematical and Physical Theories of Higher Geodesy, Vol 1, (Teubner, Leipzig, 1880), Chaps. 5–7.
- T. Vincenty,
Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations,
Survey Review
**23**(176), 88–93 (1975). - J. Danielsen,
The area under the geodesic, Survey Review
**30**(232), 61–66 (1989). - C. F. F. Karney,
Algorithms for geodesics, J. Geodesy
**87**(1) 43–55 (2013); addenda. - C. F. F. Karney, {@https://arxiv.org/abs/1102.1215v1 Geodesics on an ellipsoid of revolution}, Feb. 2011; errata.
- A geodesic bibliography.
- The wikipedia page, Geodesics on an ellipsoid.