Geodesics on an arbitrary ellipsoid of revolution

TL;DR

This work extends robust geodesic calculations to ellipsoids of revolution with arbitrary eccentricity by recasting direct/inverse problems in Legendre and Cayley elliptic integrals and by computing geodesic-polygon areas via a discrete sine transform. It introduces a stable longitude formulation, a 1D root-finding strategy for the inverse problem, and a DST-based approach for the area integral, achieving near double-precision accuracy while remaining practical for terrestrial ellipsoids. Implemented in GeographicLib with open-source availability, the methods are validated across oblate and prolate cases and tested on large polygons, including the Poland boundary, with careful attention to numerical precision and performance. The two-part methodology broadens geodesy’s applicability to highly eccentric bodies and provides a robust, efficient toolkit for precise geodesic and polygon-area computations, while clarifying when legacy Taylor-series methods remain preferable for small flattening scenarios.

Abstract

The algorithms given in Karney, J. Geodesy 87, 43-55 (2013), to compute geodesics on terrestrial ellipsoids are extended to apply to ellipsoids of revolution with arbitrary eccentricity. For the direct and inverse geodesic problems, this entails implementing the formulation in terms of elliptic integrals given by Legendre and Cayley. The integral for the area of geodesic polygons is computed in terms of the discrete sine transform of the integrand. In all cases, accuracy close to full machine precision is achieved. An open-source implementation of these algorithms is available.

Figures (6)

• Figure 1: A simple closed geodesic for an ellipsoid with $b/a=\frac{1}{4}$. The geodesic starts on the equator with azimuth $\alpha_0 \approx 51.24052^\circ$ and completes 2 complete undulations about the equator before returning to the starting longitude: (a) a top view; (b) a view at an inclination of $25^\circ$ to the equatorial plane; "N" marks the north pole. The solid and dotted lines show the visible and hidden portions of the geodesic.
• Figure 2: The solution to the hybrid problem for (a) $\phi_1 = \phi_2 = 0$ and $b/a = {\frac{1}{2}}$ (marked with crosses); (b) $\phi_1 = -30^\circ$, $\phi_2 = 20^\circ$, and $b/a = 2$ (marked with circles). This illustrates the functional dependence of $\lambda_{12}(\alpha^*_1)$.
• Figure 3: The dependence of the errors in the individual Fourier coefficients, $P_l$, on $l$ for $N=30$. The crosses, resp. squares, show the errors using the DST method, $\delta\hat{P}_l^{(N)}$, resp. the Taylor series method, $\delta\tilde{P}_l^{(N)}$. The case illustrated is for a geodesic on an ellipsoid with $b/a=\frac{1}{4}$, $\alpha_0=\frac{1}{4}\pi$. The line shows the exact value of $\lvert P_l\rvert$.
• Figure 4: The dependence of the overall error in the Fourier series approximation to $p(\sigma)$ on $N$. The crosses, resp. squares, show the errors for the DST method, $\Delta \hat{P}^{(N)}$, resp. the Taylor series method, $\Delta \tilde{P}^{(N)}$. This is for the same case as Fig. \ref{['deltaP-l']}, namely $b/a=\frac{1}{4}$, $\alpha_0=\frac{1}{4}\pi$. The horizontal line indicates the roundoff limit for double precision $2^{-53} \approx 10^{-16}$.
• Figure 5: The heavy line shows the minimum number of points $N$ required in the DST to ensure double precision accuracy as a function of $n$ for $\lvert n\rvert \le 0.99$ and for arbitrary $\alpha_0$. If we limit the possible values of $N$ to $2\times 2^j$ and $3\times 2^j$ with integer $j \ge 1$, then the result is the staircase shown as a light line.