Geodesics on an Arbitrary Ellipsoid of Revolution

Introduction

In an earlier paper, Algorithms for geodesics (Karney, 2013, henceforth referred to as AG), I presented algorithms for solving geodesic problems on an ellipsoid of revolution. This built on the classic work of Bessel (1825) and Helmert (1880) to give full double precision accuracy for ellipsoids with flattening satisfying |f| ≤ 1 100 . A small improvement, described in Appendix A.1, extended this range |f| ≤ 1
50 . Perhaps the most consequential innovation of AG was the reliable solution of the inverse geodesic problem, computing the shortest path between two arbitrary points on the ellipsoid. (Previously, the state of the art was provided by Vincenty (1975), which fails to converge in some cases.) Another important advance was turning the formulation of Danielsen (1989) for the area between a geodesic segment and the equator into an algorithm for computing accurately the area of any geodesic polygon, a polygon whose edges are geodesics. These algorithms allow geodesic problems to be reliably solved and this, in turn, has meant that they have been widely adopted in a variety of disciplines.

The present paper seeks to extend the treatment to arbitrary ellipsoid (in this paper, the term “ellipsoid” should be understood to mean “ellipsoid of revolution”). The condition |f| ≤ 1 50 covers all terrestrial applications. However, there are planetary applications where f lies outside this range (for Saturn, we have f ≈ 1 10 ). For geodesic calculations, the solutions of the direct and inverse problems, this is a matter of realizing the formulation in terms of elliptic integrals, originally given by Legendre, as a usable algorithm. Generalizing
the computation of the area of a geodesic polygon required a novel approach using the discrete sine transform. The appendices cover various small improvements in the original algorithms given in AG and compare the algorithms presented here with recent work by
Nowak and Nowak Da Costa (2022).