GeographicLib
1.40

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.
See also
Go to
NOTES
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.
Fig. 1
Fig. 1: The ellipsoidal grid. The blue (resp. green) lines are lines of constant \(\beta\) (resp. \(\omega\)); the grid spacing is 10°. 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^2b^2)\sin^2\omega+(b^2c^2)\cos^2\beta} &= \frac{b^2\sin^2\beta+c^2\cos^2\beta} {a^2b^2\sin^2\betac^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\omegac^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 nonuniform limit, with the position of the umbilic points depending on the ratio \((ab)/(bc)\).
Interconversions between the three different latitudes and longitudes and the cartesian coordinates are simple algebraic exercises.
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^2c^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^2b^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 start at an umbilical points, we have \(\gamma = 0\) and \(\delta\) determines the direction at the umbilical point. Incidentally the constant \(\gamma\) may be expressed as
\[ \gamma = (b^2c^2)\cos^2\beta\sin^2\alpha(a^2b^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}{(b^2c^2)\cos^2\beta + (a^2b^2)\sin^2\omega} &= \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^2c^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^2b^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^2c^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^2b^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^2c^2)}\). The resulting integrals 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 \(\tanh^{1}\sin\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.
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.
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\).
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.
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\).
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.
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.
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 offdiagonal 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.)
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)\).
(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\).
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).
The shortest path found by this method is unique unless:
In Vorlesungen über Dynamik, §28, Jacobi gives the following conformal mapping of the triaxial ellipsoid onto a plane
\[ \begin{align} x &= 2 \int \frac {\sqrt{a^2 \sin^2\omega + b^2 \cos^2\omega}} {\sqrt{a^2 \sin^2\omega + b^2 \cos^2\omega  c^2}}\, d\omega, \\ y &= 2 \int \frac {\sqrt{b^2 \sin^2\beta + c^2 \cos^2\beta}} {\sqrt{a^2  b^2 \sin^2\beta  c^2 \cos^2\beta}}\, d\beta. \end{align} \]
The scale of the projection is
\[ k = \frac2 {\sqrt{a^2 \sin^2\omega + b^2 (\cos^2\omega\sin^2\beta)  c^2 \cos^2\beta}}. \]
Nyrtsov, et al.,
express the projection in terms of elliptic integrals,
\[ \begin{align} x&=\frac{2b}{\sqrt{a^2c^2}} G\biggl(\omega',\frac{a^2b^2}{a^2c^2}, \frac cb\frac{\sqrt{a^2b^2}}{\sqrt{a^2c^2}}\biggr),\\ y&=\frac{2c}{\sqrt{a^2b^2}} G\biggl(\beta',\frac{b^2c^2}{a^2b^2}, i\frac ac\frac{\sqrt{b^2c^2}}{\sqrt{a^2b^2}}\biggr), \end{align} \]
where
\[ \begin{align} \tan\omega' &= \frac{\sqrt{a^2c^2}}{\sqrt{b^2c^2}} \tan\omega, \\ \tan\beta' &= \frac{\sqrt{a^2b^2}}{\sqrt{a^2c^2}} \tan\beta, \end{align} \]
and
\[ \begin{align} G(\phi,\alpha^2,k) &= \int_0^\phi \frac{\sqrt{1  k^2\sin^2\theta}}{1  \alpha^2\sin^2\theta}\,d\theta\\ &=\frac{k^2}{\alpha^2}F(\phi, k) +\biggl(1\frac{k^2}{\alpha^2}\biggr)\Pi(\phi, \alpha^2, k) \end{align} \]
is the combination of elliptic integrals of the first and third kinds that appears in the integral, given by Legendre (1811), for the longitude of a geodesic on an ellipsoid of revolution. The function \( G(\phi,\alpha^2,k) \) is provided by the EllipticFunction class. The expression for \( x \) is essentially the same at that given by Nyrtsov, et al.; I have put their expression for \( y \) into a somewhat simpler form using an imaginary modulus.
Using http://dlmf.nist.gov/19.7.E9 the projection may also be written as
\[ \begin{align} x &= \frac{2a^2}{b{\sqrt{a^2c^2}}} \Pi\biggl(\omega', \frac{a^2b^2}{b^2}, \frac cb\frac{\sqrt{a^2b^2}}{\sqrt{a^2c^2}}\biggr) \\ &\quad{}2\tan^{1}\biggl(\frac{(a^2b^2) \sin\omega' \cos\omega'} {\sqrt{a^2(b^2c^2)\sin^2\omega' + b^2(a^2c^2)\cos^2\omega'}}\biggr),\\ y &= \frac{2b^2}{c\sqrt{a^2b^2}} \Pi\biggl(\beta',\frac{b^2c^2}{c^2}, i\frac ac\frac{\sqrt{b^2c^2}}{\sqrt{a^2b^2}}\biggr) \\ &\quad{}2\sinh^{1}\biggl(\frac{(b^2c^2)\sin\beta \cos\beta} {\sqrt{b^2(a^2b^2)\sin^2\beta + c^2(a^2c^2)\cos^2\beta}}\biggr). \end{align} \]
In the case of a nearly spherical ellipsoid, the parameter of the elliptic integral (its second argument) is small.
Notes:
\[ \begin{align} x &= \frac{2a}{\sqrt{a^2  c^2}} \omega,\\ y &= \frac{2a}{\sqrt{a^2  c^2}} \int \frac {\sqrt{a^2 \sin^2\beta + c^2 \cos^2\beta}}{a\cos\beta}\, d\beta\\ &= \frac{2a}{\sqrt{a^2  c^2}} \sinh^{1}\biggl(\frac ac \tan\beta\biggr)  2\sinh^{1}\biggl(\frac{\sqrt{a^2c^2}}c\sin\beta\biggr),\\ k &= \frac{2a}{\sqrt{a^2  c^2}}\frac1{a\cos\beta}; \end{align} \]
recall that \(\beta\) is the parametric latitude in this limit.\[ f(z;e) = \cosh^{1}(z/e)  \log(2/e). \]
For \(e > 0\), this function has two square root singularities at \(\pm e\), corresponding to the two northern umbilic points. Plotting contours of its real (resp. imaginary) part gives the behavior of the lines of constant latitude (resp. longitude) near the north pole in Fig. 1. If we pass to the limit \(e\rightarrow 0\), then \( f(z;e)\rightarrow\log z\), and the two branch points merge yielding a stronger (logarithmic) singularity at \(z = 0\), concentric circles of latitude, and radial lines of longitude.\[ \begin{align} x\bigl({\textstyle\frac12}\pi\bigr) &= \frac{2a^2}{b{\sqrt{a^2c^2}}} \Pi\biggl(\frac{a^2b^2}{b^2}, \frac cb\frac{\sqrt{a^2b^2}}{\sqrt{a^2c^2}}\biggr),\\ y\bigl({\textstyle\frac12}\pi\bigr) &= \frac{2b^2}{c\sqrt{a^2b^2}} \Pi\biggl(\frac{b^2c^2}{c^2}, i\frac ac\frac{\sqrt{b^2c^2}}{\sqrt{a^2b^2}}\biggr). \end{align} \]
\[ \begin{align} a&=(6378137+35)\,\mathrm m,\\ b&=(637813735)\,\mathrm m,\\ c&=6356752\,\mathrm m,\\ \end{align} \]
we have\[ \begin{align} x\bigl({\textstyle\frac12}\pi\bigr) &= 38.396508 = 1.0000005\,q,\\ y\bigl({\textstyle\frac12}\pi\bigr) &= 103.717719 = 2.7012293\,q,\\ \end{align} \]
where\[ q = \pi \frac{\sqrt{a^2+b^2}}{\sqrt{a^2+b^22c^2}}. \]
An implementation of the Jacobi conformal projection is given here