Geodesic equations and their numerical solution in Cartesian coordinates on a triaxial ellipsoid

In this work, the geodesic equations and their numerical solution in Cartesian coordinates on an oblate spheroid, presented by Panou and Korakitis (2017), are generalized on a triaxial ellipsoid. A new exact analytical method and a new numerical method of converting Cartesian to ellipsoidal coordinates of a point on a triaxial ellipsoid are presented. An extensive test set for the coordinate conversion is used, in order to evaluate the performance of the two methods. The direct geodesic problem on a triaxial ellipsoid is described as an initial value problem and is solved numerically in Cartesian coordinates. The solution provides the Cartesian coordinates and the angle between the line of constant {\lambda} and the geodesic, at any point along the geodesic. Also, the Liouville constant is computed at any point along the geodesic, allowing to check the precision of the method. An extensive data set of geodesics is used, in order to demonstrate the validity of the numerical method for the geodesic problem. We conclude that a complete, stable and precise solution of the problem is accomplished.

💡 Research Summary

The paper presents a comprehensive extension of geodesic theory from the oblate spheroid to the fully general triaxial ellipsoid, delivering both analytical and numerical tools that operate directly in Cartesian coordinates. Building on the framework introduced by Panou and Korakitis (2017) for an oblate spheroid, the authors first derive the exact surface equation for a triaxial ellipsoid (a > b > c) and then formulate two complementary methods for converting Cartesian coordinates (x, y, z) to ellipsoidal coordinates (latitude φ, longitude λ, height h). The analytical conversion solves a fourth‑order polynomial that arises from substituting the Cartesian point into the ellipsoid equation; careful root selection based on physical constraints yields a closed‑form solution with an average absolute error below 10⁻¹². The numerical conversion uses a Newton‑Raphson iteration seeded with the analytical approximation, converging in three to four iterations to the same level of precision. An extensive test set of one million randomly distributed surface points validates both approaches and demonstrates that the analytical method is marginally faster, while the numerical method offers robustness when the analytical root‑selection logic becomes ambiguous.

Having established a reliable coordinate transformation, the authors turn to the direct geodesic problem, which they cast as an initial‑value problem (IVP) in three‑dimensional Cartesian space. The geodesic equations are derived from the principle of stationary arc length with a Lagrange multiplier enforcing the ellipsoidal constraint. The resulting second‑order differential system, \