Triaxial ellipsoids¶

Jump to:

Coordinate systems
Coordinate system conversions
- Solving for h
- Solving for u
Geodesics
Utilities
References

Introduction¶

This page provides some implementation details for the triaxial class introduced in version 2.2 (released 2024-04-09) of the Octave/MATLAB package GeographicLib (this link includes instructions on how to download and install the package). This involves a lot of new code so… Expect there to be errors in the documentation (please report). Also be prepared for the interface to change. The class includes

the solution of the direct and inverse geodesic problems,

conversions between various coordinate systems,

random sampling on the ellipsoid,

functions to aid plotting curves on the ellipsoid.

Once the package is in your path, type triaxial.doc to get an overview of the class. You can then use, e.g., help triaxial.distance, to get the documentation on specific routines. Use triaxial.demo to list various demonstrations of the geodesic capabilities.

The package works equally well in Octave and in MATLAB. However note that the solution of the geodesic problems is about 40 times faster in MATLAB. Version 2.3 (released 2024-07-09) fixes some bugs in this class.

Coordinate systems¶

A triaxial ellipsoid is the surface defined by

\[S(\mathbf r) = \frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} - 1 = 0,\]

where \(a\), \(b\), and \(c\) are the semiaxes of the ellipsoid. We will take \(a \ge b \ge c > 0\). Note that the ellipsoid is centered at the origin and that its principal axes are aligned with the coordinate axes.

A direction on the surface can be specified by unit cartesian vector tangential the surface

where \(\mathbf U(\mathbf r) = \frac12 \nabla\mathbf S(\mathbf r)\) is the normal to the surface.

A point on the surface is specified by a latitude and longitude. The geodetic latitude and longitude \((\phi, \lambda)\) are defined by

\[\hat{\mathbf U} = (\cos\phi \cos\lambda, \cos\phi \sin\lambda, \sin\phi),\]

where the “hat” symbol denote a unit vector. The parametric latitude and longitude \((\phi', \lambda')\) are defined by

\[\begin{split}x &= a \cos\phi' \cos\lambda', \\ y &= b \cos\phi' \sin\lambda', \\ z &= c \sin\phi'.\end{split}\]

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

\[\hat{\mathbf r} = ( \cos\phi'' \cos\lambda'' , \cos\phi'' \sin\lambda'' , \sin\phi'' ).\]

As with ellipsoids of revolution, the geodetic, parametric, and geocentric coordinates are closely related to one another.

Finally ellipsoid latitude and longitude \((\beta, \omega)\) are defined by

\[\begin{split}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{split}\]

A notable feature of ellipsoidal coordinates is that they are orthogonal (unlike either geodetic or parametric coordinates). The ellipsoidal azimuth is then well defined. In the limit of an oblate ellipsoid of revolution, \((\beta, \omega)\) play the roles of the parametric latitude and the longitude. For a prolate ellipsoid, these two roles are switched.

There are two useful representation of arbitrary points in three-dimensional space. There first represents positions by

\[\mathbf r = \mathbf r_0 + h \hat{\mathbf U}(\mathbf r_0),\]

where \(\mathbf r_0\) is the closest point on the ellipsoid and \(h\) is the height above the ellipsoid. Since geodetic coordinates specify the direction of \(\mathbf U\), we can also represent this point be appending the height to the geodetic coordinates to give \((\phi, \lambda, h)\).

The second uses ellipsoidal coordinates. For an arbitrary point \(\mathbf r\), we seek the value of \(u\) such that

\[\frac{x^2}{u^2 + l_a^2} + \frac{y^2}{u^2 + l_b^2} + \frac{z^2}{u^2} = 1,\]

where

\[l_a = \sqrt{a^2 - c^2}, \quad l_b = \sqrt{b^2 - c^2}\]

are linear eccentricities. This is an ellipsoid which is confocal to the original one (with semiaxes \(a, b, c\)) and whose minor semiaxis is \(u\).

Coordinate system conversions¶

Conversions between these coordinate on the surface of the ellipsoid are algebraic exercises. For example, the conversion from cartesian to geodetic coordinates proceeds as follows

\[\begin{split}\xi &= x/a^2, \quad \eta = y/b^2, \quad \zeta = z/c^2, \\ \phi &= \tan^{-1} \frac\zeta{\lVert\xi, \eta\rVert}, \\ \lambda &= \tan^{-1} \frac\eta\xi,\end{split}\]

where the quadrant of the angles should be determined by the signs of the numerators and denominators separately, using, for example, the library function atan2, and \(\lVert x, y, \ldots \rVert = \sqrt{x^2 + y^2 + \ldots}\)

Solving for \(h\)¶

Following [Ligas12], we have

\[\begin{split}\mathbf r_0 &= (x_0, y_0, z_0) = \biggl( \frac{a^2x}{p + l_a^2}, \frac{b^2y}{p + l_b^2}, \frac{c^2z}{p} \biggr), \\ h &= \hat{\mathbf U}(\mathbf r_0) \cdot (\mathbf r - \mathbf r_0),\end{split}\]

where \(p\) is the largest root of

\[f(p) = \biggl(\frac{ax}{p + l_a^2}\biggr)^2 + \biggl(\frac{by}{p + l_b^2}\biggr)^2 + \biggl(\frac{cz}{p}\biggr)^2 - 1 = 0.\]

[Ligas12] uses Newton’s method to find this root; however, with his choice of starting guess, this sometimes fails to converge. [Panou+Korakitis22] cure this defect by using the bisection method to find the root. This is guaranteed to converge but at the high computation cost of requiring many iterations.

It turns out we can easily fix the problems with Newton’s method. First of all, we note that \(f(p)\) has positive double poles at \(p = 0\), \(-l_b^2\), and \(-l_a^2\) and that \(f(p) \rightarrow -1\) for \(p \rightarrow \pm\infty\). (For now, we assume that \(x, y, z\) are all nonzero.). Therefore \(f(p)\) has a single root for \(p \in (0, \infty)\). In this region \(f'(p) < 0\) and \(f''(p) > 0\), and, as a consequence, picking a starting guess for Newton’s method between \(p = 0\) and the actual root is guaranteed to converge.

To obtain a reasonably tight bound on the root, we use

\[\begin{split}f(p) &\le \biggl(\frac{cz}{p}\biggr)^2 - 1, \\ f(p) &\le \biggl(\frac{\lVert by, cz\rVert}{p + l_b^2}\biggr)^2 - 1, \\ f(p) &\le \biggl(\frac{\lVert ax, by, cz\rVert}{p + l_a^2}\biggr)^2 - 1, \\ f(p) &\ge \biggl(\frac{\lVert ax, by, cz\rVert}{p}\biggr)^2 - 1.\end{split}\]

Because \(f'(p) < 0\) for \(p > 0\), this leads to bounds on the positive root, \(p_{\mathrm{min}} \le p \le p_{\mathrm{max}}\), where

\[\begin{split}p_{\mathrm{min}} &= \max(\lvert cz\rvert, \lVert by, cz\rVert - l_b^2, \lVert ax, by, cz\rVert - l_a^2), \\ p_{\mathrm{max}} &= \lVert ax, by, cz\rVert.\end{split}\]

[Panou+Korakitis22] substitute \(p_{\mathrm{min}} = c \lvert z\rvert\); they would get better performance using the tighter bound given here. [Ligas12] uses \(p_0 = c\lVert x, y, z\rVert\) for his initial guess; because \(f(p_0)\) can then be negative, Newton’s method may fail to converge.

In implementing Newton’s method, we neglect any term in the definition of \(f(p)\) if its numerator vanishes (even though the denominator might also vanish).

Provided that \(f(p_{\mathrm{min}}) > 0\), we can now start Newton’s method with \(p_0 = p_{\mathrm{min}}\) and this converges to the root from below. If \(f(p_{\mathrm{min}}) \le 0\) (which can only happen if \(z = 0\)), the required solution is \(p = 0\). In this case, the expression for \(\mathbf r_0\) is indeterminate, and we proceed as follows:

If \(x_0\) is indeterminate, substitute \(x_0 = 0\) (this can only happen with \(x = 0\) on a sphere).
If \(y_0\) is indeterminate, substitute \(y_0 = 0\) (this can only happen with \(y = 0\) on an oblate spheroid).
If \(z_0\) is indeterminate, substitute \(z_0 = c \sqrt{1 - x^2/a^2 - y^2/b^2}\).

A few other points to note:

This prescription obviates the need to enumerate and solve various subcases as [Panou+Korakitis22] do.
Newton’s method requires about 8 times fewer iterations compared with the bisection method.
The independent variable \(f(p)\) is shifted with respect to the one used by [Ligas12] and [Panou+Korakitis22]. This gives higher precision close to the singularity at \(p = 0\).
We accumulate the terms in \(f(p)\) in a two-word accumulator to improved the accuracy near its root.
To avoid potentially singular behavior, we initially “flush” tiny values of the components of \(\mathbf r\) to zero.
In the case where \(z_0\) is indeterminate, the sign of \(z\) should be used to determine the sign of the square root above.
If need be, this method is easily generalized to ellipsoids in higher dimensions.

Solving for \(u\)¶

Writing \(u^2 = q\), we need to find the roots of

\[g(q) = \frac{x^2}{q + l_a^2} + \frac{y^2}{q + l_b^2} + \frac{z^2}{q} - 1 = 0.\]

The structure of \(g(q)\) is very similar to \(f(p)\). Since \(g(q)\) has 3 simple poles with positive coefficients, there are three real roots and, because the rightmost pole is at \(q = 0\), just one of them is positive. As before, bounds can be put on this root \(q_{\mathrm{min}} \le q \le q_{\mathrm{max}}\), where

\[\begin{split}q_{\mathrm{min}} &= \max(z^2, y^2 + z^2 - l_b^2, x^2 + y^2 + z^2 - l_a^2), \\ q_{\mathrm{max}} &= x^2 + y^2 + z^2.\end{split}\]

As before, in implementing Newton’s method, we neglect any term in the definition of \(g(q)\) if its numerator vanishes (even though the denominator might also vanish).

Provided that \(g(q_{\mathrm{min}}) > 0\), we can now start Newton’s method with \(q_0 = q_{\mathrm{min}}\) and this converges to the root from below. If \(g(q_{\mathrm{min}}) \le 0\) (which can only happen if \(z = 0\)), the required solution is \(q = u = 0\).

Of course, we can expand out \(g(q)\) to obtain a cubic polynomial in \(q\) which cab be solved analytically, see [DLMF], Secs. 1.11(iii) and 4.43. This method is advocated by [Panou+Korakitis21]. However, this solution suffers from roundoff error when the coefficient of \(q\) is positive; in this case, the polynomial in \(1/q\) should be solved instead. Even so, the solution may be subject to unacceptable roundoff error; it may be refined by using as the starting point, \(q_0\), for Newton’s method. If \(g(q_0) < 0\), then \(q_1\) should be replaced by \(\max(q_1, q_{\mathrm{min}})\). Typically only 3–4 iterations are needed to refine the solution.

Note: tighter bounds can be placed on \(q\) using

\[\begin{split}g(q) &\le \frac{y^2}{q + l_b^2} + \frac{z^2}{q} - 1 \\ g(q) &\le \frac{x^2+y^2}{q + l_a^2} + \frac{z^2}{q} - 1 \\ g(q) &\le \frac{x^2}{q + l_a^2} + \frac{y^2+z^2}{q + l_b^2} - 1 \\ g(q) &\ge \frac{x^2+y^2}{q + l_b^2} + \frac{z^2}{q} - 1 \\ g(q) &\ge \frac{x^2}{q + l_a^2} + \frac{y^2+z^2}{q} - 1\end{split}\]

and solving the resulting quadratic equations. This yields only a marginal improvement given that we’re starting with the root of the cubic.

Geodesics¶

The problem of geodesics on a triaxial ellipsoid was solved by [Jacobi39] who reduced the problem to quadrature. Even without evaluating the integrals, this solution allowed the various properties of geodesics to be found. For an overview, see [GeographicLib-triaxial].

Explicit evaluation of Jacobi’s integrals was carried out by hand by [Cayley72] and, more recently, by [Baillard15]. Accurate evaluation of the integrals involves changing the variable of integration using elliptic integrals and elliptic functions. Unfortunately, Octave/MATLAB has poor support for these special functions, so for this implementation of the geodesic routines, I instead integrate the geodesic equations in cartesian coordinates, following [Panou+Korakitis19].

The direct problem¶

The equation for geodesics on a surface is the same as for the motion of a particle constrained to move on the surface but subject to no other forces. The centrifugal acceleration of the particle is \(-(v^2/R)\hat{\mathbf U}\) where \(R\) is the radius of curvature in the direction of \(\mathbf v\). We will take the speed to be unity (and, of course, the speed is a constant in this problem); thus time can be replaced by \(s\), the displacement along the geodesic, as the independent variable. The differential equation for the geodesic is

\[\begin{split}d\mathbf r / ds &= \mathbf v, \\ d\mathbf v / ds &= \mathbf A, \\ d^2 m / ds^2 &= -K m,\end{split}\]

where

\[\begin{split}\mathbf A &= \frac{\mathbf U}{U^2} \biggl( \frac{v_x^2}{a^2} + \frac{v_y^2}{b^2} + \frac{v_z^2}{c^2} \biggr),\\ K &= \frac1{a^2b^2c^2 U^4}.\end{split}\]

It is simplest to expression \(\mathbf r\) and \(\mathbf v\) is cartesian coordinates, since then there are no singularities in the representation.

Here \(m\) is the reduced length and \(K\) is the Gaussian curvature. It’s not necessary to determine \(m\) to solve the direct problem; however, it is an important aspect of solving the inverse problem.

We use the ODE routines provided with Octave and MATLAB to solve these equations. To make the control of the error simpler, we first scale the ellipsoid so that its middle semiaxis \(b = 1\); then all the dependent variables are of order unity. The ODE solvers take care of picking the appropriate step size for integration. In addition, they allow intermediate points on the path to be found inexpensively by polynomial interpolation.

The demonstrations triaxial.demo(n) for n = 1:5 and n = 11:15 show the result of solution of the direct problem of various initial conditions. These illustrate the distinctive properties of geodesics, i.e., that the undulate between either lines of constant \(\beta\) or lines of constant \(\omega\). In the limiting case, the geodesic repeatedly returns to opposite umbilical points.

Note well: Octave is about 40 slower than MATLAB at solving the ODEs.

The inverse problem¶

[Panou13] and [Baillard15] both attempt to solve the inverse problem, finding the shortest path between two points. However, neither offers a complete solution. A reliable method of solving the problem is obtained using the same basic method give by [Karney13] for solving the problem on an oblate ellipsoid; this is outlined in [GeographicLib-triaxial]. The key observation by [Itoh+Kiyohara04] is that the cut locus for geodesics emanating from a given point is a segment of the line of the opposite ellipsoidal latitude; see triaxial.demo(6).

The solution in the general case, involves starting with the point with the large absolute latitude, varying the azimuth at this point and find the longitude where this geodesic intersects the line of latitude for the other point. This makes use of the ability for the ODE solvers in Octave/MATLAB to stop at the occurrence of certain “events”. The azimuth can be corrected using Newton’s method (this is where the reduced length \(m\) is needed) to find the azimuth where the longitude matches that of the other point.

About 6 iterations are required for random pairs of points on a terrestrial ellipsoid. Based on the [Geodesic-testset], the overall accuracy is probably about 10 μm. The method is somewhat fragile in that it expects geodesics to behave in the way dictated by Jacobi’s solution; however, the ODE solver cannot guarantee that this is so. However by setting reasonably tight error tolerances are set on the ODE solver and deploying some other defensive tricks, the method works as long as the ellipsoid is not too eccentric. (To be safe, the ellipsoid should satisfy \(a/b \le 2\) and \(b/c \le 2\). Also avoid ellipsoids which are nearly but not quite ellipsoids of revolution; triaxial models of the earth are fine, but expect problems if the difference in the equatorial semiaxes is 1 μm.)

This method therefore provides a “working” solution of the inverse problem. A “complete” solution will involve using Jacobi’s solution. This will remove the sloppiness involved in using an ODE solver. An initial implementation of Jacobi’s solution was used to create the [Geodesic-testset].

Jacobi’s solution¶

I have coded up Jacobi’s solution to the direct problem in MATLAB using the [Chebfun] package. This allows the indefinite integrals in Jacobi’s solution to be evaluated accurately. I do not include this functionality in the triaxial class because

Chebfun is not compatible with Octave;
MATLAB’s support for elliptic integrals and elliptic functions with modulus close to 1 is deficient — this leads to inaccuracies for geodesics which graze the umbilical points.

I will reimplement the solution in the C++ version of GeographicLib. This will make more consistent use of Fourier series (in contrast, Chebfun switches to a Chebyshev series when asked to integrate a Fourier series) and use GeographicLib’s implementation of elliptic integrals and elliptic functions.

With this in place, the solution of the inverse problem should be straightforward. Jacobi does not include an expression for the reduced length \(m\), so I will use some method other than Newton’s for finding the azimuth, such as the method of [Chandrupatla97].

Utilities¶

You can sample points (and directions) uniformly on the ellipsoid with cart2rand, see [Marples+Williams23]

The function horizon returns points on the horizon of the ellipsoid when viewed from a distant viewpoint in the direction \(\mathbf V\). These points satisfy

\[\begin{split}\mathbf U \cdot \mathbf V &= 0\\ \biggl(\frac x{a^2}, \frac y{b^2}, \frac z{c^2}\biggl) \cdot \mathbf V&= 0\\ \biggl(\frac xa, \frac yb, \frac zc\biggl) \cdot \biggl(\frac{V_x}a, \frac{V_y}b, \frac{V_z}c\biggl) &= 0\end{split}\]

The first vector in the last equation gives points on a unit sphere, and these are on the horizon of the sphere when viewed from the direction given by the second vector. So the ellipsoidal horizon is obtained by computing this spherical horizon (a circle) and scaling the cartesian components by \((a, b, c)\).