Geodesics on a triaxial ellipsoid (original) (raw)

Back to Finding nearest neighbors. Forward to Jacobi's conformal projection. Up to Contents.

Jacobi (1839) showed that the problem of geodesics on a triaxial ellipsoid (with 3 unequal axes) can be reduced to quadrature. Despite this, the detailed behavior of the geodesics is not very well known. In this section, I briefly give Jacobi's solution and illustrate the behavior of the geodesics and outline an algorithm for the solution of the inverse problem.

Triaxial coordinate systems

Consider the ellipsoid defined by

\[ f = \frac{X^2}{a^2} + \frac{Y^2}{b^2} + \frac{Z^2}{c^2} = 1, \]

where, without loss of generality, \( a \ge b \ge c \gt 0\). A point on the surface is specified by a latitude and longitude. The geographical latitude and longitude \((\phi, \lambda)\) are defined by

\[ \frac{\nabla f}{\left| \nabla f\right|} = \left( \begin{array}{c} \cos\phi \cos\lambda \\ \cos\phi \sin\lambda \\ \sin\phi \end{array}\right). \]

The parametric latitude and longitude \((\phi', \lambda')\) are defined by

\[ \begin{align} X &= a \cos\phi' \cos\lambda', \\ Y &= b \cos\phi' \sin\lambda', \\ Z &= c \sin\phi'. \end{align} \]

Jacobi employed the ellipsoidal latitude and longitude \((\beta, \omega)\) defined by

\[ \begin{align} X &= a \cos\omega \frac{\sqrt{a^2 - b^2\sin^2\beta - c^2\cos^2\beta}} {\sqrt{a^2 - c^2}}, \\ Y &= b \cos\beta \sin\omega, \\ Z &= c \sin\beta \frac{\sqrt{a^2\sin^2\omega + b^2\cos^2\omega - c^2}} {\sqrt{a^2 - c^2}}. \end{align} \]

Grid lines of constant \(\beta\) and \(\omega\) are given in Fig. 1.

Geodesic grid on a triaxial ellipsoid Fig. 1

Fig. 1: The ellipsoidal grid. The blue (resp. green) lines are lines of constant \(\beta\) (resp. \(\omega\)); the grid spacing is 15°. Also shown in red are two of the principal sections of the ellipsoid, defined by \(X = 0\) and \(Z = 0\). The third principal section, \(Y = 0\), is covered by the lines \(\beta = \pm 90^\circ\) and \(\omega = 90^\circ \pm 90^\circ\). These lines meet at four umbilical points (two of which are visible in this figure) where the principal radii of curvature are equal. The parameters of the ellipsoid are \(a = 1.01\), \(b = 1\), \(c = 0.8\), and it is viewed in an orthographic projection from a point above \(\phi = 40^\circ\), \(\lambda = 30^\circ\). These parameters were chosen to accentuate the ellipsoidal effects on geodesics (relative to those on the earth) while still allowing the connection to an ellipsoid of revolution to be made.

The grid lines of the ellipsoid coordinates are "lines of curvature" on the ellipsoid, i.e., they are parallel to the direction of principal curvature (Monge, 1796). They are also intersections of the ellipsoid with confocal systems of hyperboloids of one and two sheets (Dupin, 1813). Finally they are geodesic ellipses and hyperbolas defined using two adjacent umbilical points. For example, the lines of constant \(\beta\) in Fig. 1 can be generated with the familiar string construction for ellipses with the ends of the string pinned to the two umbilical points.

The element of length on the ellipsoid in ellipsoidal coordinates is given by

\[ \begin{align} \frac{ds^2}{(a^2-b^2)\sin^2\omega+(b^2-c^2)\cos^2\beta} &= \frac{b^2\sin^2\beta+c^2\cos^2\beta} {a^2-b^2\sin^2\beta-c^2\cos^2\beta} d\beta^2 \\ &\qquad+ \frac{a^2\sin^2\omega+b^2\cos^2\omega} {a^2\sin^2\omega+b^2\cos^2\omega-c^2} d\omega^2. \end{align} \]

The torus \((\omega, \beta) \in [-\pi,\pi] \times [-\pi,\pi]\) covers the ellipsoid twice. In order to facilitate passing to the limit of an oblate ellipsoid, we may regard as the principal sheet \([-\pi,\pi] \times [-\frac12\pi,\frac12\pi]\) and insert branch cuts at \(\beta=\pm\frac12\pi\). The rule for switching sheets is

\[ \begin{align} \omega & \rightarrow -\omega,\\ \beta & \rightarrow \pi-\beta,\\ \alpha & \rightarrow \pi+\alpha, \end{align} \]

where \(\alpha\) is the heading of a path, relative to a line of constant \(\omega\).

In the limit \(b\rightarrow a\) (resp. \(b\rightarrow c\)), the umbilic points converge on the \(Z\) (resp. \(X\)) axis and an oblate (resp. prolate) ellipsoid is obtained with \(\beta\) (resp. \(\omega\)) becoming the standard parametric latitude and \(\omega\) (resp. \(\beta\)) becoming the standard longitude. The sphere is a non-uniform limit, with the position of the umbilic points depending on the ratio \((a-b)/(b-c)\).

Inter-conversions between the three different latitudes and longitudes and the cartesian coordinates are simple algebraic exercises.

Jacobi's solution

Solving the geodesic problem for an ellipsoid of revolution is, from the mathematical point of view, trivial; because of symmetry, geodesics have a constant of the motion (analogous to the angular momentum) which was found by Clairaut (1733). By 1806 (with the work of Legendre, Oriani, et al.), there was a complete understanding of the qualitative behavior of geodesics on an ellipsoid of revolution.

On the other hand, geodesics on a triaxial ellipsoid have no obvious constant of the motion and thus represented a challenging "unsolved" problem in the first half of the nineteenth century. Jacobi discovered that the geodesic equations are separable if they are expressed in ellipsoidal coordinates. You can get an idea of the importance Jacobi attached to his discovery from the letter he wrote to his friend and neighbor Bessel:

The day before yesterday, I reduced to quadrature the problem of geodesic lines on an ellipsoid with three unequal axes. They are the simplest formulas in the world, Abelian integrals, which become the well known elliptic integrals if 2 axes are set equal.
Königsberg, 28th Dec. '38.

On the same day he wrote a similar letter to the editor of Compte Rendus and his result was published in J. Crelle in (1839) with a French translation (from German) appearing in J. Liouville in (1841).

Here is the solution, exactly as given by Jacobi here (with minor changes in notation):

\[ \begin{align} \delta &= \int \frac {\sqrt{b^2\sin^2\beta + c^2\cos^2\beta}\,d\beta} {\sqrt{a^2 - b^2\sin^2\beta - c^2\cos^2\beta} \sqrt{(b^2-c^2)\cos^2\beta - \gamma}}\\ &\quad - \int \frac {\sqrt{a^2\sin^2\omega + b^2\cos^2\omega}\,d\omega} {\sqrt{a^2\sin^2\omega + b^2\cos^2\omega - c^2} \sqrt{(a^2-b^2)\sin^2\omega + \gamma}} \end{align} \]

As Jacobi notes "a function of the angle \(\beta\) equals a function of the angle \(\omega\). These two functions are just Abelian integrals…" Two constants \(\delta\) and \(\gamma\) appear in the solution. Typically \(\delta\) is zero if the lower limits of the integrals are taken to be the starting point of the geodesic and the direction of the geodesics is determined by \(\gamma\). However for geodesics that go through an umbilical point, we have \(\gamma = 0\) and \(\delta\) determines the direction at the umbilical point. Incidentally the constant \(\gamma\) may be expressed as

\[ \gamma = (b^2-c^2)\cos^2\beta\sin^2\alpha-(a^2-b^2)\sin^2\omega\cos^2\alpha \]

where \(\alpha\) is the angle the geodesic makes with lines of constant \(\omega\). In the limit \(b\rightarrow a\), this reduces to \(\cos\beta\sin\alpha = \text{const.}\), the familiar Clairaut relation. A nice derivation of Jacobi's result is given by Darboux (1894) §§583–584 where he gives the solution found by Liouville (1846) for general quadratic surfaces. In this formulation, the distance along the geodesic, \(s\), is also found using

\[ \begin{align} \frac{ds}{(a^2-b^2)\sin^2\omega + (b^2-c^2)\cos^2\beta} &= \frac {\sqrt{b^2\sin^2\beta + c^2\cos^2\beta}\,d\beta} {\sqrt{a^2 - b^2\sin^2\beta - c^2\cos^2\beta} \sqrt{(b^2-c^2)\cos^2\beta - \gamma}}\\ &= \frac {\sqrt{a^2\sin^2\omega + b^2\cos^2\omega}\,d\omega} {\sqrt{a^2\sin^2\omega + b^2\cos^2\omega - c^2} \sqrt{(a^2-b^2)\sin^2\omega + \gamma}} \end{align} \]

An alternative expression for the distance is

\[ \begin{align} ds &= \frac {\sqrt{b^2\sin^2\beta + c^2\cos^2\beta} \sqrt{(b^2-c^2)\cos^2\beta - \gamma}\,d\beta} {\sqrt{a^2 - b^2\sin^2\beta - c^2\cos^2\beta}}\\ &\quad {}+ \frac {\sqrt{a^2\sin^2\omega + b^2\cos^2\omega} \sqrt{(a^2-b^2)\sin^2\omega + \gamma}\,d\omega} {\sqrt{a^2\sin^2\omega + b^2\cos^2\omega - c^2}} \end{align} \]

Jacobi's solution is a convenient way to compute geodesics on an ellipsoid. Care must be taken with the signs of the square roots (which are determined by the initial azimuth of the geodesic). Also if \(\gamma \gt 0\) (resp. \(\gamma \lt 0\)), then the \(\beta\) (resp. \(\omega\)) integrand diverges. The integrand can be transformed into a finite one by a change of variable, e.g., \(\sin\beta = \sin\sigma \sqrt{1 - \gamma/(b^2-c^2)}\). The resulting integrands are periodic, so the behavior of an arbitrarily long geodesic is entirely captured by tabulating the integrals over a single period.

The situation is more complicated if \(\gamma = 0\) (corresponding to umbilical geodesics). Both integrands have simple poles at the umbilical points. However, this behavior may be subtracted from the integrands to yield (for example) the sum of a term involving \(\sinh^{-1}\tan\beta\) and a finite integral. Since both integrals contain similar logarithmic singularities they can be equated (thus fixing the ratio \(\cos\beta/\sin\omega\) at the umbilical point) and connection formulas can be found which allow the geodesic to be followed through the umbilical point. The study of umbilical geodesics was of special interest to a group of Irish mathematicians in the 1840's and 1850's, including Michael and William Roberts (twins!), Hart, Graves, and Salmon.

Survey of triaxial geodesics

Before delving into the nature of geodesics on a triaxial geodesic, it is worth reviewing geodesics on an ellipsoid of revolution. There are two classes of simple closed geodesics (i.e., geodesics which close on themselves without intersection): the equator and all the meridians. All other geodesics oscillate between two equal and opposite circles of latitude; but after completing a full oscillation in latitude these fall slightly short (for an oblate ellipsoid) of completing a full circuit in longitude.

Turning to the triaxial case, we find that there are only 3 simple closed geodesics, the three principal sections of the ellipsoid given by \(X = 0\), \(Y = 0\), and \(Z = 0\). To survey the other geodesics, it is convenient to consider geodesics which intersect the middle principal section, \(Y = 0\), at right angles. Such geodesics are shown in Figs. 2–6, where I use the same ellipsoid parameters as in Fig. 1 and the same viewing direction. In addition, the three principal ellipses are shown in red in each of these figures.

If the starting point is \(\beta_1 \in (-90^\circ, 90^\circ)\), \(\omega_1 = 0\), and \(\alpha_1 = 90^\circ\), then the geodesic encircles the ellipsoid in a "circumpolar" sense. The geodesic oscillates north and south of the equator; on each oscillation it completes slightly less that a full circuit around the ellipsoid resulting in the geodesic filling the area bounded by the two latitude lines \(\beta = \pm \beta_1\). Two examples are given in Figs. 2 and 3. Figure 2 shows practically the same behavior as for an oblate ellipsoid of revolution (because \(a \approx b\)). However, if the starting point is at a higher latitude (Fig. 3) the distortions resulting from \(a \ne b\) are evident.

Example of a circumpolar geodesic on a
triaxial ellipsoid Fig. 2

Fig. 2: Example of a circumpolar geodesic on a triaxial ellipsoid. The starting point of this geodesic is \(\beta_1 = 45.1^\circ\), \(\omega_1 = 0^\circ\), and \(\alpha_1 = 90^\circ\).

Another example of a circumpolar geodesic on a
triaxial ellipsoid Fig. 3

Fig. 3: Another example of a circumpolar geodesic on a triaxial ellipsoid. The starting point of this geodesic is \(\beta_1 = 87.48^\circ\), \(\omega_1 = 0^\circ\), and \(\alpha_1 = 90^\circ\).

If the starting point is \(\beta_1 = 90^\circ\), \(\omega_1 \in (0^\circ, 180^\circ)\), and \(\alpha_1 = 180^\circ\), then the geodesic encircles the ellipsoid in a "transpolar" sense. The geodesic oscillates east and west of the ellipse \(X = 0\); on each oscillation it completes slightly more that a full circuit around the ellipsoid resulting in the geodesic filling the area bounded by the two longitude lines \(\omega = \omega_1\) and \(\omega = 180^\circ - \omega_1\). If \(a = b\), all meridians are geodesics; the effect of \(a \ne b\) causes such geodesics to oscillate east and west. Two examples are given in Figs. 4 and 5.

Example of a transpolar geodesic on a
triaxial ellipsoid Fig. 4

Fig. 4: Example of a transpolar geodesic on a triaxial ellipsoid. The starting point of this geodesic is \(\beta_1 = 90^\circ\), \(\omega_1 = 39.9^\circ\), and \(\alpha_1 = 180^\circ\).

Another example of a transpolar geodesic on a
triaxial ellipsoid Fig. 5

Fig. 5: Another example of a transpolar geodesic on a triaxial ellipsoid. The starting point of this geodesic is \(\beta_1 = 90^\circ\), \(\omega_1 = 9.966^\circ\), and \(\alpha_1 = 180^\circ\).

If the starting point is \(\beta_1 = 90^\circ\), \(\omega_1 = 0^\circ\) (an umbilical point), and \(\alpha_1 = 135^\circ\) (the geodesic leaves the ellipse \(Y = 0\) at right angles), then the geodesic repeatedly intersects the opposite umbilical point and returns to its starting point. However on each circuit the angle at which it intersects \(Y = 0\) becomes closer to \(0^\circ\) or \(180^\circ\) so that asymptotically the geodesic lies on the ellipse \(Y = 0\). This is shown in Fig. 6. Note that a single geodesic does not fill an area on the ellipsoid.

Example of an umbilical geodesic on a
triaxial ellipsoid Fig. 6

Fig. 6: Example of an umbilical geodesic on a triaxial ellipsoid. The starting point of this geodesic is \(\beta_1 = 90^\circ\), \(\omega_1 = 0^\circ\), and \(\alpha_1 = 135^\circ\) and the geodesics is followed forwards and backwards until it lies close to the plane \(Y = 0\) in both directions.

Umbilical geodesics enjoy several interesting properties.

Through any point on the ellipsoid, there are two umbilical geodesics.
The geodesic distance between opposite umbilical points is the same regardless of the initial direction of the geodesic.
Whereas the closed geodesics on the ellipses \(X = 0\) and \(Z = 0\) are stable (an geodesic initially close to and nearly parallel to the ellipse remains close to the ellipse), the closed geodesic on the ellipse \(Y = 0\), which goes through all 4 umbilical points, is unstable. If it is perturbed, it will swing out of the plane \(Y = 0\) and flip around before returning to close to the plane. (This behavior may repeat depending on the nature of the initial perturbation.).

The stability of closed geodesics

The stability of the three simple closed geodesics can be determined by examining the properties of Jacobi's solution. In particular the unstable behavior of umbilical geodesics was shown by Hart (1849). However an alternative approach is to use the equations that Gauss (1828) gives for a perturbed geodesic

\[ \frac {d^2m}{ds^2} + Km = 0 \]

where \(m\) is the distance of perturbed geodesic from a reference geodesic and \(K\) is the Gaussian curvature of the surface. If the reference geodesic is closed, then this is a linear homogeneous differential equation with periodic coefficients. In fact it's a special case of Hill's equation which can be treated using Floquet theory, see DLMF, §28.29. Using the notation of §3 of Algorithms for geodesics, the stability is determined by computing the reduced length \(m_{12}\) and the geodesic scales \(M_{12}, M_{21}\) over half the perimeter of the ellipse and determining the eigenvalues \(\lambda_{1,2}\) of

\[ {\cal M} = \left(\begin{array}{cc} M_{12} & m_{12}\\ -\frac{1 - M_{12}M_{21}}{m_{12}} & M_{21} \end{array}\right). \]

Because \(\mathrm{det}\,{\cal M} = 1\), the eigenvalues are determined by \(\mathrm{tr}\,{\cal M}\). In particular if \(\left|\mathrm{tr}\,{\cal M}\right| < 2\), we have \(\left|\lambda_{1,2}\right| = 1\) and the solution is stable; if \(\left|\mathrm{tr}\,{\cal M}\right| > 2\), one of \(\left|\lambda_{1,2}\right|\) is larger than unity and the solution is (exponentially) unstable. In the transition case, \(\left|\mathrm{tr}\,{\cal M}\right| = 2\), the solution is stable provided that the off-diagonal elements of \({\cal M}\) are zero; otherwise the solution is linearly unstable.

The exponential instability of the geodesic on the ellipse \(Y = 0\) is confirmed by this analysis and results from the resonance between the natural frequency of the equation for \(m\) and the driving frequency when \(b\) lies in \((c, a)\). If \(b\) is equal to either of the other axes (and the triaxial ellipsoid degenerates to an ellipsoid of revolution), then the solution is linearly unstable. (For example, a geodesic is which is close to a meridian on an oblate ellipsoid, slowly moves away from that meridian.)

The inverse problem

In order to solve the inverse geodesic problem, it helps to have an understanding of the properties of all the geodesics emanating from a single point \((\beta_1, \omega_1)\).

If the point is an umbilical point, all the lines meet at the opposite umbilical point.
Otherwise, the first envelope of the geodesics is a 4-pointed astroid. The cusps of the astroid lie on either \(\beta = - \beta_1\) or \(\omega = \omega_1 + \pi\); see Sinclair (2003).
All geodesics intersect (or, in the case of \(\alpha_1 = 0\) or \(\pi\), touch) the line \(\omega = \omega_1 + \pi\).
All geodesics intersect (or, in the case of \(\alpha_1 = \pm\pi/2\), touch) the line \(\beta = -\beta_1\).
Two geodesics with azimuths \(\pm\alpha_1\) first intersect on \(\omega = \omega_1 + \pi\) and their lengths to the point of intersection are equal.
Two geodesics with azimuths \(\alpha_1\) and \(\pi-\alpha_1\) first intersect on \(\beta = -\beta_1\) and their lengths to the point of intersection are equal.

(These assertions follow directly from the equations for the geodesics; some of them are somewhat surprising given the asymmetries of the ellipsoid.) Consider now terminating the geodesics from \((\beta_1, \omega_1)\) at the point where they first intersect (or touch) the line \(\beta = -\beta_1\). To focus the discussion, take \(\beta_1 \le 0\).

The geodesics completely fill the portion of the ellipsoid satisfying \(\beta \le -\beta_1\).
None of geodesics intersect any other geodesics.
Any initial portion of these geodesics is a shortest path.
Each geodesic intersects the line \(\beta = \beta_2\), where \(\beta_1 < \beta_2 < -\beta_1\), exactly once.
For a given \(\beta_2\), this defines a continuous monotonic mapping of the circle of azimuths \(\alpha_1\) to the circle of longitudes \(\omega_2\).
If \(\beta_2 = \pm \beta_1\), then the previous two assertions need to be modified similarly to the case for an ellipsoid of revolution.

These properties show that the inverse problem can be solved using techniques similar to those employed for an ellipsoid of revolution (see §4 of Algorithms for geodesics).

If the points are opposite umbilical points, an arbitrary \(\alpha_1\) may be chosen.
If the points are neighboring umbilical points, the shortest path lies on the ellipse \(Y = 0\).
If only one point is an umbilicial point, the azimuth at the non-umbilical point is found using the generalization of Clairaut's equation (given above) with \(\gamma = 0\).
Treat the cases where the geodesic might follow a line of constant \(\beta\). There are two such cases: (a) the points lie on the ellipse \(Z = 0\), \(\beta = 0\), on a general ellipsoid and (b) the points lie on an ellipse whose major axis is the \(X\) axis, \(\beta\) arbitrary, on a prolate ellipsoid ( \(a > b = c\)). Determine the reduced length \(m_{12}\) for the geodesic which is the shorter path along the ellipse. If \(m_{12} \ge 0\), then this is the shortest path on the ellipsoid; otherwise proceed to the general case (next).
Swap the points, if necessary, so that the first point is the one closest to a pole. Estimate \(\alpha_1\) (by some means) and solve the hybrid problem, i.e., determine the longitude \(\omega_2\) corresponding to the first intersection of the geodesic with \(\beta = \beta_2\). Adjust \(\alpha_1\) so that the value of \(\omega_2\) matches the given \(\omega_2\) (there is a single root). If a sufficiently close solution can be found, Newton's method can be employed since the necessary derivative can be expressed in terms of the reduced length \(m_{12}\).

The shortest path found by this method is unique unless:

The length of the geodesic vanishes \(s_{12}=0\), in which case any constant can be added to the azimuths.
The points are opposite umbilical points. In this case, \(\alpha_1\) can take on any value and \(\alpha_2\) needs to be adjusted to maintain the value of \(\tan\alpha_1 / \tan\alpha_2\). Note that \(\alpha\) increases by \(\pm 90^\circ\) as the geodesic passes through an umbilical point, depending on whether the geodesic is considered as passing to the right or left of the point. Here \(\alpha_2\) is the forward azimuth at the second umbilical point, i.e., its azimuth immediately after passage through the umbilical point.
\(\beta_1 + \beta_2 = 0\) and \(\cos\alpha_1\) and \(\cos\alpha_2\) have opposite signs. In this case, there another shortest geodesic with azimuths \(\pi - \alpha_1\) and \(\pi - \alpha_2\).

I have implemented this solution to the inverse geodesic problem on a triaxial ellipsoid in Octave/MATLAB. However, this is uses the method of Panou and Korakitis (2019) to solve the direct problem while searching for the inverse solution.

Jacobi's conformal projection

This material is now on its own page; see Jacobi's conformal projection.

Back to Finding nearest neighbors. Forward to Jacobi's conformal projection. Up to Contents.