CINXE.COM

<!DOCTYPE html> <html lang="en" data-content_root="../"> <head> <meta charset="utf-8" /> <meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="viewport" content="width=device-width, initial-scale=1" /> <title>Triaxial ellipsoids — GeographicLib 2.3 documentation</title> <link rel="stylesheet" type="text/css" href="../_static/pygments.css?v=fa44fd50" /> <link rel="stylesheet" type="text/css" href="../_static/classic.css?v=def86cc0" /> <script src="../_static/documentation_options.js?v=57236720"></script> <script src="../_static/doctools.js?v=9a2dae69"></script> <script src="../_static/sphinx_highlight.js?v=dc90522c"></script> <script async="async" src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script> <link rel="index" title="Index" href="../genindex.html" /> <link rel="search" title="Search" href="../search.html" /> <link rel="prev" title="Research" href="research.html" /> </head><body> <div class="related" role="navigation" aria-label="Related"> <h3>Navigation</h3> <ul> <li class="right" style="margin-right: 10px"> <a href="../genindex.html" title="General Index" accesskey="I">index</a></li> <li class="right" > <a href="research.html" title="Research" accesskey="P">previous</a> |</li> <li class="nav-item nav-item-0"><a href="../index.html">GeographicLib 2.3 documentation</a> »</li> <li class="nav-item nav-item-this"><a href="">Triaxial ellipsoids</a></li> </ul> </div> <div class="document"> <div class="documentwrapper"> <div class="bodywrapper"> <div class="body" role="main"> <section id="triaxial-ellipsoids"> <span id="triaxial"></span><h1>Triaxial ellipsoids<a class="headerlink" href="#triaxial-ellipsoids" title="Link to this heading">¶</a></h1> <p>Jump to:</p> <ul class="simple"> <li><p><a class="reference internal" href="#tricoords"><span class="std std-ref">Coordinate systems</span></a></p></li> <li><dl class="simple"> <dt><a class="reference internal" href="#tricoordconv"><span class="std std-ref">Coordinate system conversions</span></a></dt><dd><ul> <li><p><a class="reference internal" href="#tricoordconvh"><span class="std std-ref">Solving for h</span></a></p></li> <li><p><a class="reference internal" href="#tricoordconvu"><span class="std std-ref">Solving for u</span></a></p></li> </ul> </dd> </dl> </li> <li><dl class="simple"> <dt><a class="reference internal" href="#trigeod"><span class="std std-ref">Geodesics</span></a></dt><dd><ul> <li><p><a class="reference internal" href="#trigeoddirect"><span class="std std-ref">The direct problem</span></a></p></li> <li><p><a class="reference internal" href="#trigeodinverse"><span class="std std-ref">The inverse problem</span></a></p></li> <li><p><a class="reference internal" href="#trigeodjac"><span class="std std-ref">Jacobi’s solution</span></a></p></li> </ul> </dd> </dl> </li> <li><p><a class="reference internal" href="#triutils"><span class="std std-ref">Utilities</span></a></p></li> <li><p><a class="reference internal" href="#trirefs"><span class="std std-ref">References</span></a></p></li> </ul> <section id="introduction"> <span id="triintro"></span><h2>Introduction<a class="headerlink" href="#introduction" title="Link to this heading">¶</a></h2> <p>This page provides some implementation details for the <code class="docutils literal notranslate"><span class="pre">triaxial</span></code> class introduced in version 2.2 (released 2024-04-09) of the Octave/MATLAB package <a class="reference external" href="https://github.com/geographiclib/geographiclib-octave#readme">GeographicLib</a> (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</p> <blockquote> <div><ul class="simple"> <li><p>the solution of the direct and inverse geodesic problems,</p></li> <li><p>conversions between various coordinate systems,</p></li> <li><p>random sampling on the ellipsoid,</p></li> <li><p>functions to aid plotting curves on the ellipsoid.</p></li> </ul> </div></blockquote> <p>Once the package is in your path, type <code class="docutils literal notranslate"><span class="pre">triaxial.doc</span></code> to get an overview of the class. You can then use, e.g., <code class="docutils literal notranslate"><span class="pre">help</span> <span class="pre">triaxial.distance</span></code>, to get the documentation on specific routines. Use <code class="docutils literal notranslate"><span class="pre">triaxial.demo</span></code> to list various demonstrations of the geodesic capabilities.</p> <p>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.</p> </section> <section id="coordinate-systems"> <span id="tricoords"></span><h2>Coordinate systems<a class="headerlink" href="#coordinate-systems" title="Link to this heading">¶</a></h2> <p>A triaxial ellipsoid is the surface defined by</p> <div class="math notranslate nohighlight"> \[S(\mathbf r) = \frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} - 1 = 0,\]</div> <p>where <span class="math notranslate nohighlight">\(a\)</span>, <span class="math notranslate nohighlight">\(b\)</span>, and <span class="math notranslate nohighlight">\(c\)</span> are the semiaxes of the ellipsoid. We will take <span class="math notranslate nohighlight">\(a \ge b \ge c > 0\)</span>. Note that the ellipsoid is centered at the origin and that its principal axes are aligned with the coordinate axes.</p> <p>A direction on the surface can be specified by unit cartesian vector tangential the surface</p> <p>where <span class="math notranslate nohighlight">\(\mathbf U(\mathbf r) = \frac12 \nabla\mathbf S(\mathbf r)\)</span> is the normal to the surface.</p> <p>A point on the surface is specified by a latitude and longitude. The <em>geodetic</em> latitude and longitude <span class="math notranslate nohighlight">\((\phi, \lambda)\)</span> are defined by</p> <div class="math notranslate nohighlight"> \[\hat{\mathbf U} = (\cos\phi \cos\lambda, \cos\phi \sin\lambda, \sin\phi),\]</div> <p>where the “hat” symbol denote a unit vector. The <em>parametric</em> latitude and longitude <span class="math notranslate nohighlight">\((\phi', \lambda')\)</span> are defined by</p> <div class="math notranslate nohighlight"> \[\begin{split}x &= a \cos\phi' \cos\lambda', \\ y &= b \cos\phi' \sin\lambda', \\ z &= c \sin\phi'.\end{split}\]</div> <p>The <em>geocentric</em> latitude and longitude <span class="math notranslate nohighlight">\((\phi'', \lambda'')\)</span> are defined by</p> <div class="math notranslate nohighlight"> \[\hat{\mathbf r} = ( \cos\phi'' \cos\lambda'' , \cos\phi'' \sin\lambda'' , \sin\phi'' ).\]</div> <p>As with ellipsoids of revolution, the geodetic, parametric, and geocentric coordinates are closely related to one another.</p> <p>Finally <em>ellipsoid</em> latitude and longitude <span class="math notranslate nohighlight">\((\beta, \omega)\)</span> are defined by</p> <div class="math notranslate nohighlight"> \[\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}\]</div> <p>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, <span class="math notranslate nohighlight">\((\beta, \omega)\)</span> play the roles of the parametric latitude and the longitude. For a prolate ellipsoid, these two roles are switched.</p> <p>There are two useful representation of arbitrary points in three-dimensional space. There first represents positions by</p> <div class="math notranslate nohighlight"> \[\mathbf r = \mathbf r_0 + h \hat{\mathbf U}(\mathbf r_0),\]</div> <p>where <span class="math notranslate nohighlight">\(\mathbf r_0\)</span> is the closest point on the ellipsoid and <span class="math notranslate nohighlight">\(h\)</span> is the height above the ellipsoid. Since geodetic coordinates specify the direction of <span class="math notranslate nohighlight">\(\mathbf U\)</span>, we can also represent this point be appending the height to the geodetic coordinates to give <span class="math notranslate nohighlight">\((\phi, \lambda, h)\)</span>.</p> <p>The second uses ellipsoidal coordinates. For an arbitrary point <span class="math notranslate nohighlight">\(\mathbf r\)</span>, we seek the value of <span class="math notranslate nohighlight">\(u\)</span> such that</p> <div class="math notranslate nohighlight"> \[\frac{x^2}{u^2 + l_a^2} + \frac{y^2}{u^2 + l_b^2} + \frac{z^2}{u^2} = 1,\]</div> <p>where</p> <div class="math notranslate nohighlight"> \[l_a = \sqrt{a^2 - c^2}, \quad l_b = \sqrt{b^2 - c^2}\]</div> <p>are linear eccentricities. This is an ellipsoid which is confocal to the original one (with semiaxes <span class="math notranslate nohighlight">\(a, b, c\)</span>) and whose minor semiaxis is <span class="math notranslate nohighlight">\(u\)</span>.</p> <section id="coordinate-system-conversions"> <span id="tricoordconv"></span><h3>Coordinate system conversions<a class="headerlink" href="#coordinate-system-conversions" title="Link to this heading">¶</a></h3> <p>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</p> <div class="math notranslate nohighlight"> \[\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}\]</div> <p>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 <span class="math notranslate nohighlight">\(\lVert x, y, \ldots \rVert = \sqrt{x^2 + y^2 + \ldots}\)</span></p> </section> <section id="solving-for-h"> <span id="tricoordconvh"></span><h3>Solving for <span class="math notranslate nohighlight">\(h\)</span><a class="headerlink" href="#solving-for-h" title="Link to this heading">¶</a></h3> <p>Following <a class="reference internal" href="#ligas12" id="id1"><span>[Ligas12]</span></a>, we have</p> <div class="math notranslate nohighlight"> \[\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}\]</div> <p>where <span class="math notranslate nohighlight">\(p\)</span> is the largest root of</p> <div class="math notranslate nohighlight"> \[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.\]</div> <p><a class="reference internal" href="#ligas12" id="id2"><span>[Ligas12]</span></a> uses Newton’s method to find this root; however, with his choice of starting guess, this sometimes fails to converge. <a class="reference internal" href="#panou-korakitis22" id="id3"><span>[Panou+Korakitis22]</span></a> 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.</p> <p>It turns out we can easily fix the problems with Newton’s method. First of all, we note that <span class="math notranslate nohighlight">\(f(p)\)</span> has positive double poles at <span class="math notranslate nohighlight">\(p = 0\)</span>, <span class="math notranslate nohighlight">\(-l_b^2\)</span>, and <span class="math notranslate nohighlight">\(-l_a^2\)</span> and that <span class="math notranslate nohighlight">\(f(p) \rightarrow -1\)</span> for <span class="math notranslate nohighlight">\(p \rightarrow \pm\infty\)</span>. (For now, we assume that <span class="math notranslate nohighlight">\(x, y, z\)</span> are all nonzero.). Therefore <span class="math notranslate nohighlight">\(f(p)\)</span> has a single root for <span class="math notranslate nohighlight">\(p \in (0, \infty)\)</span>. In this region <span class="math notranslate nohighlight">\(f'(p) < 0\)</span> and <span class="math notranslate nohighlight">\(f''(p) > 0\)</span>, and, as a consequence, picking a starting guess for Newton’s method between <span class="math notranslate nohighlight">\(p = 0\)</span> and the actual root is guaranteed to converge.</p> <p>To obtain a reasonably tight bound on the root, we use</p> <div class="math notranslate nohighlight"> \[\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}\]</div> <p>Because <span class="math notranslate nohighlight">\(f'(p) < 0\)</span> for <span class="math notranslate nohighlight">\(p > 0\)</span>, this leads to bounds on the positive root, <span class="math notranslate nohighlight">\(p_{\mathrm{min}} \le p \le p_{\mathrm{max}}\)</span>, where</p> <div class="math notranslate nohighlight"> \[\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}\]</div> <p><a class="reference internal" href="#panou-korakitis22" id="id4"><span>[Panou+Korakitis22]</span></a> substitute <span class="math notranslate nohighlight">\(p_{\mathrm{min}} = c \lvert z\rvert\)</span>; they would get better performance using the tighter bound given here. <a class="reference internal" href="#ligas12" id="id5"><span>[Ligas12]</span></a> uses <span class="math notranslate nohighlight">\(p_0 = c\lVert x, y, z\rVert\)</span> for his initial guess; because <span class="math notranslate nohighlight">\(f(p_0)\)</span> can then be negative, Newton’s method may fail to converge.</p> <p>In implementing Newton’s method, we neglect any term in the definition of <span class="math notranslate nohighlight">\(f(p)\)</span> if its numerator vanishes (even though the denominator might also vanish).</p> <p>Provided that <span class="math notranslate nohighlight">\(f(p_{\mathrm{min}}) > 0\)</span>, we can now start Newton’s method with <span class="math notranslate nohighlight">\(p_0 = p_{\mathrm{min}}\)</span> and this converges to the root from below. If <span class="math notranslate nohighlight">\(f(p_{\mathrm{min}}) \le 0\)</span> (which can only happen if <span class="math notranslate nohighlight">\(z = 0\)</span>), the required solution is <span class="math notranslate nohighlight">\(p = 0\)</span>. In this case, the expression for <span class="math notranslate nohighlight">\(\mathbf r_0\)</span> is indeterminate, and we proceed as follows:</p> <ul class="simple"> <li><p>If <span class="math notranslate nohighlight">\(x_0\)</span> is indeterminate, substitute <span class="math notranslate nohighlight">\(x_0 = 0\)</span> (this can only happen with <span class="math notranslate nohighlight">\(x = 0\)</span> on a sphere).</p></li> <li><p>If <span class="math notranslate nohighlight">\(y_0\)</span> is indeterminate, substitute <span class="math notranslate nohighlight">\(y_0 = 0\)</span> (this can only happen with <span class="math notranslate nohighlight">\(y = 0\)</span> on an oblate spheroid).</p></li> <li><p>If <span class="math notranslate nohighlight">\(z_0\)</span> is indeterminate, substitute <span class="math notranslate nohighlight">\(z_0 = c \sqrt{1 - x^2/a^2 - y^2/b^2}\)</span>.</p></li> </ul> <p>A few other points to note:</p> <ul class="simple"> <li><p>This prescription obviates the need to enumerate and solve various subcases as <a class="reference internal" href="#panou-korakitis22" id="id6"><span>[Panou+Korakitis22]</span></a> do.</p></li> <li><p>Newton’s method requires about 8 times fewer iterations compared with the bisection method.</p></li> <li><p>The independent variable <span class="math notranslate nohighlight">\(f(p)\)</span> is shifted with respect to the one used by <a class="reference internal" href="#ligas12" id="id7"><span>[Ligas12]</span></a> and <a class="reference internal" href="#panou-korakitis22" id="id8"><span>[Panou+Korakitis22]</span></a>. This gives higher precision close to the singularity at <span class="math notranslate nohighlight">\(p = 0\)</span>.</p></li> <li><p>We accumulate the terms in <span class="math notranslate nohighlight">\(f(p)\)</span> in a two-word accumulator to improved the accuracy near its root.</p></li> <li><p>To avoid potentially singular behavior, we initially “flush” tiny values of the components of <span class="math notranslate nohighlight">\(\mathbf r\)</span> to zero.</p></li> <li><p>In the case where <span class="math notranslate nohighlight">\(z_0\)</span> is indeterminate, the sign of <span class="math notranslate nohighlight">\(z\)</span> should be used to determine the sign of the square root above.</p></li> <li><p>If need be, this method is easily generalized to ellipsoids in higher dimensions.</p></li> </ul> </section> <section id="solving-for-u"> <span id="tricoordconvu"></span><h3>Solving for <span class="math notranslate nohighlight">\(u\)</span><a class="headerlink" href="#solving-for-u" title="Link to this heading">¶</a></h3> <p>Writing <span class="math notranslate nohighlight">\(u^2 = q\)</span>, we need to find the roots of</p> <div class="math notranslate nohighlight"> \[g(q) = \frac{x^2}{q + l_a^2} + \frac{y^2}{q + l_b^2} + \frac{z^2}{q} - 1 = 0.\]</div> <p>The structure of <span class="math notranslate nohighlight">\(g(q)\)</span> is very similar to <span class="math notranslate nohighlight">\(f(p)\)</span>. Since <span class="math notranslate nohighlight">\(g(q)\)</span> has 3 simple poles with positive coefficients, there are three real roots and, because the rightmost pole is at <span class="math notranslate nohighlight">\(q = 0\)</span>, just one of them is positive. As before, bounds can be put on this root <span class="math notranslate nohighlight">\(q_{\mathrm{min}} \le q \le q_{\mathrm{max}}\)</span>, where</p> <div class="math notranslate nohighlight"> \[\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}\]</div> <p>As before, in implementing Newton’s method, we neglect any term in the definition of <span class="math notranslate nohighlight">\(g(q)\)</span> if its numerator vanishes (even though the denominator might also vanish).</p> <p>Provided that <span class="math notranslate nohighlight">\(g(q_{\mathrm{min}}) > 0\)</span>, we can now start Newton’s method with <span class="math notranslate nohighlight">\(q_0 = q_{\mathrm{min}}\)</span> and this converges to the root from below. If <span class="math notranslate nohighlight">\(g(q_{\mathrm{min}}) \le 0\)</span> (which can only happen if <span class="math notranslate nohighlight">\(z = 0\)</span>), the required solution is <span class="math notranslate nohighlight">\(q = u = 0\)</span>.</p> <p>Of course, we can expand out <span class="math notranslate nohighlight">\(g(q)\)</span> to obtain a cubic polynomial in <span class="math notranslate nohighlight">\(q\)</span> which cab be solved analytically, see <a class="reference internal" href="#dlmf" id="id9"><span>[DLMF]</span></a>, Secs. 1.11(iii) and 4.43. This method is advocated by <a class="reference internal" href="#panou-korakitis21" id="id10"><span>[Panou+Korakitis21]</span></a>. However, this solution suffers from roundoff error when the coefficient of <span class="math notranslate nohighlight">\(q\)</span> is positive; in this case, the polynomial in <span class="math notranslate nohighlight">\(1/q\)</span> 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, <span class="math notranslate nohighlight">\(q_0\)</span>, for Newton’s method. If <span class="math notranslate nohighlight">\(g(q_0) < 0\)</span>, then <span class="math notranslate nohighlight">\(q_1\)</span> should be replaced by <span class="math notranslate nohighlight">\(\max(q_1, q_{\mathrm{min}})\)</span>. Typically only 3–4 iterations are needed to refine the solution.</p> <p>Note: tighter bounds can be placed on <span class="math notranslate nohighlight">\(q\)</span> using</p> <div class="math notranslate nohighlight"> \[\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}\]</div> <p>and solving the resulting quadratic equations. This yields only a marginal improvement given that we’re starting with the root of the cubic.</p> </section> </section> <section id="geodesics"> <span id="trigeod"></span><h2>Geodesics<a class="headerlink" href="#geodesics" title="Link to this heading">¶</a></h2> <p>The problem of geodesics on a triaxial ellipsoid was solved by <a class="reference internal" href="#jacobi39" id="id11"><span>[Jacobi39]</span></a> 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 <a class="reference internal" href="#geographiclib-triaxial" id="id12"><span>[GeographicLib-triaxial]</span></a>.</p> <p>Explicit evaluation of Jacobi’s integrals was carried out by hand by <a class="reference internal" href="#cayley72" id="id13"><span>[Cayley72]</span></a> and, more recently, by <a class="reference internal" href="#baillard15" id="id14"><span>[Baillard15]</span></a>. 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 <a class="reference internal" href="#panou-korakitis19" id="id15"><span>[Panou+Korakitis19]</span></a>.</p> <section id="the-direct-problem"> <span id="trigeoddirect"></span><h3>The direct problem<a class="headerlink" href="#the-direct-problem" title="Link to this heading">¶</a></h3> <p>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 <span class="math notranslate nohighlight">\(-(v^2/R)\hat{\mathbf U}\)</span> where <span class="math notranslate nohighlight">\(R\)</span> is the radius of curvature in the direction of <span class="math notranslate nohighlight">\(\mathbf v\)</span>. 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 <span class="math notranslate nohighlight">\(s\)</span>, the displacement along the geodesic, as the independent variable. The differential equation for the geodesic is</p> <div class="math notranslate nohighlight"> \[\begin{split}d\mathbf r / ds &= \mathbf v, \\ d\mathbf v / ds &= \mathbf A, \\ d^2 m / ds^2 &= -K m,\end{split}\]</div> <p>where</p> <div class="math notranslate nohighlight"> \[\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}\]</div> <p>It is simplest to expression <span class="math notranslate nohighlight">\(\mathbf r\)</span> and <span class="math notranslate nohighlight">\(\mathbf v\)</span> is cartesian coordinates, since then there are no singularities in the representation.</p> <p>Here <span class="math notranslate nohighlight">\(m\)</span> is the reduced length and <span class="math notranslate nohighlight">\(K\)</span> is the Gaussian curvature. It’s not necessary to determine <span class="math notranslate nohighlight">\(m\)</span> to solve the direct problem; however, it is an important aspect of solving the inverse problem.</p> <p>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 <span class="math notranslate nohighlight">\(b = 1\)</span>; 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.</p> <p>The demonstrations <code class="docutils literal notranslate"><span class="pre">triaxial.demo(n)</span></code> for <code class="docutils literal notranslate"><span class="pre">n</span> <span class="pre">=</span> <span class="pre">1:5</span></code> and <code class="docutils literal notranslate"><span class="pre">n</span> <span class="pre">=</span> <span class="pre">11:15</span></code> 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 <span class="math notranslate nohighlight">\(\beta\)</span> or lines of constant <span class="math notranslate nohighlight">\(\omega\)</span>. In the limiting case, the geodesic repeatedly returns to opposite umbilical points.</p> <p>Note well: Octave is about 40 slower than MATLAB at solving the ODEs.</p> </section> <section id="the-inverse-problem"> <span id="trigeodinverse"></span><h3>The inverse problem<a class="headerlink" href="#the-inverse-problem" title="Link to this heading">¶</a></h3> <p><a class="reference internal" href="#panou13" id="id16"><span>[Panou13]</span></a> and <a class="reference internal" href="#baillard15" id="id17"><span>[Baillard15]</span></a> 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 <a class="reference internal" href="#karney13" id="id18"><span>[Karney13]</span></a> for solving the problem on an oblate ellipsoid; this is outlined in <a class="reference internal" href="#geographiclib-triaxial" id="id19"><span>[GeographicLib-triaxial]</span></a>. The key observation by <a class="reference internal" href="#itoh-kiyohara04" id="id20"><span>[Itoh+Kiyohara04]</span></a> is that the <em>cut locus</em> for geodesics emanating from a given point is a segment of the line of the opposite ellipsoidal latitude; see <code class="docutils literal notranslate"><span class="pre">triaxial.demo(6)</span></code>.</p> <p>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 <span class="math notranslate nohighlight">\(m\)</span> is needed) to find the azimuth where the longitude matches that of the other point.</p> <p>About 6 iterations are required for random pairs of points on a terrestrial ellipsoid. Based on the <a class="reference internal" href="#geodesic-testset" id="id21"><span>[Geodesic-testset]</span></a>, 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 <span class="math notranslate nohighlight">\(a/b \le 2\)</span> and <span class="math notranslate nohighlight">\(b/c \le 2\)</span>. 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.)</p> <p>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 <a class="reference internal" href="#geodesic-testset" id="id22"><span>[Geodesic-testset]</span></a>.</p> </section> <section id="jacobi-s-solution"> <span id="trigeodjac"></span><h3>Jacobi’s solution<a class="headerlink" href="#jacobi-s-solution" title="Link to this heading">¶</a></h3> <p>I have coded up Jacobi’s solution to the direct problem in MATLAB using the <a class="reference internal" href="#chebfun" id="id23"><span>[Chebfun]</span></a> package. This allows the indefinite integrals in Jacobi’s solution to be evaluated accurately. I do not include this functionality in the <code class="docutils literal notranslate"><span class="pre">triaxial</span></code> class because</p> <ul class="simple"> <li><p>Chebfun is not compatible with Octave;</p></li> <li><p>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.</p></li> </ul> <p>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.</p> <p>With this in place, the solution of the inverse problem should be straightforward. Jacobi does not include an expression for the reduced length <span class="math notranslate nohighlight">\(m\)</span>, so I will use some method other than Newton’s for finding the azimuth, such as the method of <a class="reference internal" href="#chandrupatla97" id="id24"><span>[Chandrupatla97]</span></a>.</p> </section> </section> <section id="utilities"> <span id="triutils"></span><h2>Utilities<a class="headerlink" href="#utilities" title="Link to this heading">¶</a></h2> <p>You can sample points (and directions) uniformly on the ellipsoid with <code class="docutils literal notranslate"><span class="pre">cart2rand</span></code>, see <a class="reference internal" href="#marples-williams23" id="id25"><span>[Marples+Williams23]</span></a></p> <p>The function <code class="docutils literal notranslate"><span class="pre">horizon</span></code> returns points on the horizon of the ellipsoid when viewed from a distant viewpoint in the direction <span class="math notranslate nohighlight">\(\mathbf V\)</span>. These points satisfy</p> <div class="math notranslate nohighlight"> \[\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}\]</div> <p>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 <span class="math notranslate nohighlight">\((a, b, c)\)</span>.</p> </section> <section id="references"> <span id="trirefs"></span><h2>References<a class="headerlink" href="#references" title="Link to this heading">¶</a></h2> <div role="list" class="citation-list"> <div class="citation" id="baillard15" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span>Baillard15<span class="fn-bracket">]</span></span> <span class="backrefs">(<a role="doc-backlink" href="#id14">1</a>,<a role="doc-backlink" href="#id17">2</a>)</span> <p>Baillard. <a class="reference external" href="https://hp41programs.yolasite.com/geod3axial.php">Geodesics on a triaxial ellipsoid for the HP-41</a> (2015).</p> </div> <div class="citation" id="cayley72" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span><a role="doc-backlink" href="#id13">Cayley72</a><span class="fn-bracket">]</span></span> <p>Cayley, <a class="reference external" href="https://books.google.com/books?id=S4znAAAAMAAJ&pg=PA31">On the geodesic lines on an ellipsoid</a> (1872).</p> </div> <div class="citation" id="chandrupatla97" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span><a role="doc-backlink" href="#id24">Chandrupatla97</a><span class="fn-bracket">]</span></span> <p>Chandrupatla, <a class="reference external" href="https://doi.org/10.1016/s0965-9978(96)00051-8">A new hybrid quadratic/bisection algorithm for finding the zero of a nonlinear function without using derivatives</a> (1997).</p> </div> <div class="citation" id="chebfun" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span><a role="doc-backlink" href="#id23">Chebfun</a><span class="fn-bracket">]</span></span> <p>Chebfun, <a class="reference external" href="https://www.chebfun.org">Numerical computing with functions</a> (2014).</p> </div> <div class="citation" id="dlmf" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span><a role="doc-backlink" href="#id9">DLMF</a><span class="fn-bracket">]</span></span> <p>Olver et al., <a class="reference external" href="https://dlmf.nist.gov">NIST Handbook of Mathematical Functions</a> (2010).</p> </div> <div class="citation" id="geodesic-testset" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span>Geodesic-testset<span class="fn-bracket">]</span></span> <span class="backrefs">(<a role="doc-backlink" href="#id21">1</a>,<a role="doc-backlink" href="#id22">2</a>)</span> <p>Karney, <a class="reference external" href="https://doi.org/10.5281/zenodo.12510796">Test set of geodesics on a trixial ellipsoid</a> (2024).</p> </div> <div class="citation" id="geographiclib-triaxial" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span>GeographicLib-triaxial<span class="fn-bracket">]</span></span> <span class="backrefs">(<a role="doc-backlink" href="#id12">1</a>,<a role="doc-backlink" href="#id19">2</a>)</span> <p>Karney, <a class="reference external" href="https://geographiclib.sourceforge.io/1.29/triaxial.html">Geodesics on a triaxial ellipsoid</a> (2013).</p> </div> <div class="citation" id="itoh-kiyohara04" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span><a role="doc-backlink" href="#id20">Itoh+Kiyohara04</a><span class="fn-bracket">]</span></span> <p>Itoh & Kiyohara, <a class="reference external" href="https://doi.org/10.1007/s00229-004-0455-z">The cut loci and the conjugate loci on ellipsoids</a> (2004).</p> </div> <div class="citation" id="jacobi39" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span><a role="doc-backlink" href="#id11">Jacobi39</a><span class="fn-bracket">]</span></span> <p>Jacobi, <a class="reference external" href="https://doi.org/10.1515/crll.1839.19.309">Note von der geodätischen Linie auf einem Ellipsoid und den verschiedenen Anwendungen einer merkwürdigen analytischen Substitution</a> (1839).</p> </div> <div class="citation" id="karney13" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span><a role="doc-backlink" href="#id18">Karney13</a><span class="fn-bracket">]</span></span> <p>Karney, <a class="reference external" href="http://dx.doi.org/10.1007/s00190-012-0578-z">Algorithms for geodesics</a> (2013).</p> </div> <div class="citation" id="ligas12" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span>Ligas12<span class="fn-bracket">]</span></span> <span class="backrefs">(<a role="doc-backlink" href="#id1">1</a>,<a role="doc-backlink" href="#id2">2</a>,<a role="doc-backlink" href="#id5">3</a>,<a role="doc-backlink" href="#id7">4</a>)</span> <p>Ligas, <a class="reference external" href="http://dx.doi.org/10.1007/s11200-011-9017-5">Two modified algorithms to transform Cartesian to geodetic coordinates on a triaxial ellipsoid</a> (2012).</p> </div> <div class="citation" id="marples-williams23" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span><a role="doc-backlink" href="#id25">Marples+Williams23</a><span class="fn-bracket">]</span></span> <p>Marples & Williams, <a class="reference external" href="https://doi.org/10.1007/s11075-023-01628-4">Patch area and uniform sampling on the surface of any ellipsoid</a> (2023).</p> </div> <div class="citation" id="panou13" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span><a role="doc-backlink" href="#id16">Panou13</a><span class="fn-bracket">]</span></span> <p>Panou, <a class="reference external" href="10.https://doi.org/2478/jogs-2013-0028">The geodesic boundary value problem and its solution on a triaxial ellipsoid</a> (2013).</p> </div> <div class="citation" id="panou-korakitis19" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span><a role="doc-backlink" href="#id15">Panou+Korakitis19</a><span class="fn-bracket">]</span></span> <p>Panou & Korakitis, <a class="reference external" href="http://dx.doi.org/10.1515/jogs-2019-0001">Geodesic equations and their numerical solution in Cartesian coordinates on a triaxial ellipsoid</a> (2019).</p> </div> <div class="citation" id="panou-korakitis21" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span><a role="doc-backlink" href="#id10">Panou+Korakitis21</a><span class="fn-bracket">]</span></span> <p>Panou & Korakitis, <a class="reference external" href="http://dx.doi.org/10.1515/jogs-2020-0126">Analytical and numerical methods of converting Cartesian to ellipsoidal coordinates</a> (2021).</p> </div> <div class="citation" id="panou-korakitis22" role="doc-biblioentry"> <span class="label"><span class="fn-bracket">[</span>Panou+Korakitis22<span class="fn-bracket">]</span></span> <span class="backrefs">(<a role="doc-backlink" href="#id3">1</a>,<a role="doc-backlink" href="#id4">2</a>,<a role="doc-backlink" href="#id6">3</a>,<a role="doc-backlink" href="#id8">4</a>)</span> <p>Panou & Korakitis, <a class="reference external" href="http://dx.doi.org/10.1007/s00190-022-01650-9">Cartesian to geodetic coordinates conversion on a triaxial ellipsoid using the bisection method</a> (2022).</p> </div> </div> </section> </section> <div class="clearer"></div> </div> </div> </div> <div class="sphinxsidebar" role="navigation" aria-label="Main"> <div class="sphinxsidebarwrapper"> <div> <h3><a href="../index.html">Table of Contents</a></h3> <ul> <li><a class="reference internal" href="#">Triaxial ellipsoids</a><ul> <li><a class="reference internal" href="#introduction">Introduction</a></li> <li><a class="reference internal" href="#coordinate-systems">Coordinate systems</a><ul> <li><a class="reference internal" href="#coordinate-system-conversions">Coordinate system conversions</a></li> <li><a class="reference internal" href="#solving-for-h">Solving for <span class="math notranslate nohighlight">\(h\)</span></a></li> <li><a class="reference internal" href="#solving-for-u">Solving for <span class="math notranslate nohighlight">\(u\)</span></a></li> </ul> </li> <li><a class="reference internal" href="#geodesics">Geodesics</a><ul> <li><a class="reference internal" href="#the-direct-problem">The direct problem</a></li> <li><a class="reference internal" href="#the-inverse-problem">The inverse problem</a></li> <li><a class="reference internal" href="#jacobi-s-solution">Jacobi’s solution</a></li> </ul> </li> <li><a class="reference internal" href="#utilities">Utilities</a></li> <li><a class="reference internal" href="#references">References</a></li> </ul> </li> </ul> </div> <div> <h4>Previous topic</h4> <p class="topless"><a href="research.html" title="previous chapter">Research</a></p> </div> <div role="note" aria-label="source link"> <h3>This Page</h3> <ul class="this-page-menu"> <li><a href="../_sources/doc/triaxial.rst.txt" rel="nofollow">Show Source</a></li> </ul> </div> <search id="searchbox" style="display: none" role="search"> <h3 id="searchlabel">Quick search</h3> <div class="searchformwrapper"> <form class="search" action="../search.html" method="get"> <input type="text" name="q" aria-labelledby="searchlabel" autocomplete="off" autocorrect="off" autocapitalize="off" spellcheck="false"/> <input type="submit" value="Go" /> </form> </div> </search> <script>document.getElementById('searchbox').style.display = "block"</script> </div> </div> <div class="clearer"></div> </div> <div class="related" role="navigation" aria-label="Related"> <h3>Navigation</h3> <ul> <li class="right" style="margin-right: 10px"> <a href="../genindex.html" title="General Index" >index</a></li> <li class="right" > <a href="research.html" title="Research" >previous</a> |</li> <li class="nav-item nav-item-0"><a href="../index.html">GeographicLib 2.3 documentation</a> »</li> <li class="nav-item nav-item-this"><a href="">Triaxial ellipsoids</a></li> </ul> </div> <div class="footer" role="contentinfo"> © Copyright 2009–2024, Charles Karney. Created using <a href="https://www.sphinx-doc.org/">Sphinx</a> 8.0.2. </div> </body> </html>