Jacobi's solution for geodesics on a triaxial ellipsoid

TL;DR

This work delivers a robust numerical implementation of Jacobi's classical solution for geodesics on a triaxial ellipsoid, reducing the geodesic problem to the evaluation of a set of one-dimensional integrals via ellipsoidal coordinates. It employs Fourier-series representations and Clenshaw summation to achieve high-precision indefinite integrals and uses Newton-based solvers for the direct problem and Chandrupatla's method for the inverse problem, with careful handling of umbilics and special cases. The main contributions include a scalable, double-precision-capable algorithm that matches or exceeds ODE-based approaches in accuracy while offering consistent performance across geodesic lengths, and a concrete, open implementation in GeographicLib v2.7 tested on Cayley’s ellipsoid and other triaxial models. The results demonstrate high reliability and efficiency for both direct and inverse geodesic problems, making Jacobi’s solution practically viable for geodesy and planetary science applications.

Abstract

On Boxing Day, 1838, Jacobi found a solution to the problem of geodesics on a triaxial ellipsoid, with the course of the geodesic and the distance along it given in terms of one-dimensional integrals. Here, a numerical implementation of this solution is described. This entails accurately evaluating the integrals and solving the resulting coupled system of equations. The inverse problem, finding the shortest path between two points on the ellipsoid, can then be solved using a similar method as for biaxial ellipsoids.

Figures (9)

• Figure 1: The ellipsoidal grid showing lines of constant $\beta$ and $\omega$. The grid spacing is $15^\circ$. The heavy lines show the minor ($X = 0$ or $\cos\omega = 0$), median ($Y = 0$ or $\cos\beta \sin\omega = 0$), and major ($Z = 0$ or $\sin\beta = 0$) principal ellipses of the ellipsoid. The parameters of the ellipsoid are $a = 1.01$, $b = 1$, and $c = 0.8$, and it is viewed in an orthographic projection looking at the point with geodetic coordinates $\phi = 40^\circ$, $\lambda = 30^\circ$.
• Figure 2: Samples of geodesics on an ellipsoid with the same parameters and viewpoint as Fig. \ref{['graticule']}. The parameters of the geodesics are given in Table \ref{['sample-tab']}. These figures are adapted from figures that the author contributed to Wikipedia wikigeod.
• Figure 3: The graphical method for plotting umbilical geodesics given by cayley72. The semiaxes of the ellipsoid are $a = \sqrt{1000}$, $b = \sqrt{500}$, $c = \sqrt{250}$, and it is viewed here in an orthogonal projection along the $Y$ axis. The dotted lines are lines of $\beta = \beta^{(f)}_n$ (going bottom to top) or $\omega = \omega^{(f)}_n$ (going left to right), as defined in the text, with the separation constant given by $\Delta^{(f)} = 1/\sqrt{160}$ (matching Cayley's choice). In this figure, the lines of constant $\beta$ and $\omega$ are labeled with the corresponding values of $n$ (these values are given in Table \ref{['cayley-tab']}). The geodesics, shown as heavy lines, connect the vertices of the resulting mesh, and all converge on the umbilic labeled $U$, or on the neighboring one (not shown in the figure).
• Figure 4: Marking the distance along geodesics. The heavy lines are those geodesics in Fig. \ref{['cayley-fig']} that converge on the umbilic $U$. As in that figure, the dotted lines are lines of constant $\beta$ or $\omega$; however, in this case, the values are given by equal increments of $\Delta^{(g)} = 1/10$ in the distance functions (these values are given in Table \ref{['cayley-tab']}). The dashed lines connect the vertices of the mesh, and these mark off distance intervals of $b/10$ along the geodesics.
• Figure 5: Geodesics emanating from a single point on an ellipsoid. The ellipsoid and viewpoint are the same as Fig. \ref{['graticule']}. The geodesics start at $\beta_1 = -34.46^\circ$, $\omega_1 = -149.94^\circ$, and the azimuths $\alpha_1$ are multiples of $7.5^\circ$. The geodetic coordinates of the starting point are $\phi_1 = -40^\circ$, $\lambda_1 = 30^\circ-180^\circ$, the opposite of the viewing direction. Part (a) shows the geodesics followed up to the points where they are no longer the shortest geodesics. The union of such points, the cut locus, is shown as a heavy line and is a segment of the line of curvature $\beta = -\beta_1$. Part (b) shows the geodesics from (a) as dotted lines, and these are continued (as solid lines) until they meet at the line of curvature $\omega = \omega_1 + 180^\circ$ (shown as a heavy line).