[Paper Review] Geodesics on an ellipsoid of revolution

Charles F. F. Karney|arXiv (Cornell University)|Mar 1, 2009

Historical Geography and Cartography25 references21 citations

TL;DR

This paper presents highly accurate, non-iterative algorithms for solving the direct and inverse geodesic problems on an ellipsoid of revolution, achieving sub-nanometer precision. It introduces a novel formulation using the auxiliary sphere and Jacobian elliptic functions, enabling robust solutions for triangulation, map projections, maritime boundaries, and polygonal area computation with guaranteed convergence and high numerical stability.

ABSTRACT

Algorithms for the computation of the forward and inverse geodesic problems for an ellipsoid of revolution are derived. These are accurate to better than 15 nm when applied to the terrestrial ellipsoids. The solutions of other problems involving geodesics (triangulation, projections, maritime boundaries, and polygonal areas) are investigated.

Motivation & Objective

• To develop a numerically stable and highly accurate solution for the direct and inverse geodesic problems on an ellipsoid of revolution.
• To overcome the limitations of existing algorithms—particularly Vincenty’s method—by eliminating fixed-order flattening dependence and fixing convergence failure in nearly antipodal cases.
• To generalize geodesic solutions to prolate ellipsoids and extend applications to triangulation, map projections, maritime boundaries, and spherical polygon area computation.
• To provide a comprehensive theoretical and computational framework that integrates reduced length and other advanced geodetic concepts for use in GeographicLib.

Proposed method

• Formulates geodesics using the auxiliary sphere method, transforming the ellipsoidal problem into a solvable spherical triangle problem via conformal mapping.
• Employs Jacobian elliptic functions (sn, dn, am) and their integrals (F, E, G) to express geodesic arc length, azimuth, and position with high precision.
• Solves the direct problem non-iteratively using series expansions in terms of flattening $f$, third flattening $n$, or eccentricity $e^2$, ensuring accuracy to better than 15 nm.
• Applies Newton’s method with carefully chosen initial guesses to solve the inverse problem, guaranteeing convergence in a few iterations even for nearly antipodal points.
• Introduces the concept of reduced length $m_{12}$ as a key invariant for geodesic triangles, enabling accurate triangulation and area computation.
• Adapts the gnomonic projection and spherical polygon area formulas to the ellipsoidal case using geodesic distances and azimuths, preserving accuracy.

Experimental results

Research questions

• RQ1 How can the direct and inverse geodesic problems on an ellipsoid of revolution be solved with sub-nanometer accuracy and guaranteed convergence?
• RQ2 What modifications are required to extend geodesic algorithms to prolate ellipsoids, and how do the solutions differ from those for oblate ellipsoids?
• RQ3 How can geodesic-based methods be applied to compute areas of spherical polygons and define maritime boundaries with high precision?
• RQ4 What role do elliptic integrals and Jacobian functions play in enabling a non-iterative solution to the direct geodesic problem?
• RQ5 How can the reduced length concept improve the accuracy and stability of triangulation and geodetic network computations?

Key findings

• The direct geodesic problem is solved non-iteratively with accuracy better than 15 nm using series expansions in $f$, $n$, or $e^2$, valid for both oblate and prolate ellipsoids.
• The inverse geodesic problem converges reliably in a few iterations using Newton’s method, even for nearly antipodal points, overcoming a key flaw in Vincenty’s algorithm.
• The reduced length $m_{12}$ is a fundamental invariant that enables accurate triangulation and area computation on geodesic networks.
• For prolate ellipsoids, the longitude difference $\lambda_{12}$ is no longer monotonic with respect to initial azimuth $\alpha_1$, requiring modified starting point selection in Newton’s method.
• The area of a geodesic polygon is computed as $2\pi n - \sum E_{i-1,i}$, where $n$ is the number of full eastward encirclements, ensuring high-precision area integration.
• The geodesic area integral is reformulated using $G(\phi, \alpha^2, k)$ and complex elliptic functions for prolate cases, with $c^2$ defined via $\tan^{-1}(\sqrt{-e^2})/\sqrt{-e^2}$ to maintain real-valued expressions.