geographiclib-doc 1.45-2 / usr / share / doc / geographiclib / html / triaxial.html

This file is indexed.

/usr/share/doc/geographiclib/html/triaxial.html is in geographiclib-doc 1.45-2.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below. 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<meta name="generator" content="Doxygen 1.8.9.1"/>
<title>GeographicLib: Geodesics on a triaxial ellipsoid</title>
<link href="tabs.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="dynsections.js"></script>
<script type="text/x-mathjax-config">
  MathJax.Hub.Config({
    extensions: ["tex2jax.js"],
    jax: ["input/TeX","output/HTML-CSS"],
});
</script><script src="/usr/share/javascript/mathjax/MathJax.js/MathJax.js"></script>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
</head>
<body>
<div id="top"><!-- do not remove this div, it is closed by doxygen! -->
<div id="titlearea">
<table cellspacing="0" cellpadding="0">
 <tbody>
 <tr style="height: 56px;">
  <td style="padding-left: 0.5em;">
   <div id="projectname">GeographicLib
   &#160;<span id="projectnumber">1.45</span>
   </div>
  </td>
 </tr>
 </tbody>
</table>
</div>
<!-- end header part -->
<!-- Generated by Doxygen 1.8.9.1 -->
  <div id="navrow1" class="tabs">
    <ul class="tablist">
      <li><a href="index.html"><span>Main&#160;Page</span></a></li>
      <li class="current"><a href="pages.html"><span>Related&#160;Pages</span></a></li>
      <li><a href="namespaces.html"><span>Namespaces</span></a></li>
      <li><a href="annotated.html"><span>Classes</span></a></li>
      <li><a href="files.html"><span>Files</span></a></li>
    </ul>
  </div>
</div><!-- top -->
<div class="header">
  <div class="headertitle">
<div class="title">Geodesics on a triaxial ellipsoid </div>  </div>
</div><!--header-->
<div class="contents">
<div class="textblock"><center> Back to <a class="el" href="geodesic.html">Geodesics on an ellipsoid of revolution</a>. Forward to <a class="el" href="jacobi.html">Jacobi's conformal projection</a>. Up to <a class="el" href="index.html#contents">Contents</a>. </center><p>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.</p>
<p>See also</p><ul>
<li>The wikipedia page, <a href="https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid#Geodesics_on_a_triaxial_ellipsoid">Geodesics on a triaxial ellipsoid</a>.</li>
</ul>
<p>Go to</p><ul>
<li><a class="el" href="triaxial.html#triaxial-coords">Triaxial coordinate systems</a></li>
<li><a class="el" href="triaxial.html#triaxial-jacobi">Jacobi's solution</a></li>
<li><a class="el" href="triaxial.html#triaxial-survey">Survey of triaxial geodesics</a></li>
<li><a class="el" href="triaxial.html#triaxial-stab">The stability of closed geodesics</a></li>
<li><a class="el" href="triaxial.html#triaxial-inverse">The inverse problem</a></li>
</ul>
<p><b>NOTES</b></p><ol type="1">
<li>A triaxial ellipsoid approximates the earth only slightly better than an ellipsoid of revolution. If you are really considering measuring distances on the earth using a triaxial ellipsoid, you should also be worrying about the shape of the geoid, which essentially makes the geodesic problem a hopeless mess; see, for example, <a href="http://arxiv.org/abs/1112.3231">Waters (2011)</a>.</li>
<li>There is nothing new in this section. It is just an exercise in exploring Jacobi's solution. My interest here is in generating long geodesics with the correct long-time behavior. Arnold gives a nice qualitative description of the solution in <em>Mathematical Methods of Classical Mechanics</em> (2nd edition, Springer, 1989), pp. 264&ndash;266.</li>
<li>Possible reasons this problem might, nevertheless, be of interest are:<ul>
<li>It is the first example of a dynamical system which has a non-trivial constant of motion. As such, Jacobi's paper generated a lot of excitement and was followed by many papers elaborating his solution. In particular, the unstable behavior of one of the closed geodesics of the ellipsoid, is an early example of a system with a positive Lyapunov exponent (one of the essential ingredients for chaotic behavior in dynamical systems).</li>
<li>Knowledge of ellipsoidal coordinates (used by Jacobi) might be useful in other areas of geodesy.</li>
<li>Geodesics which pass through the pole on an ellipsoid of revolution represent a degenerate class (they are all closed and all pass through the opposite pole). It is of interest to see how this degeneracy is broken with a surface with a more general shape.</li>
<li>Similarly, it is of interest to see how the Mercator projection of the ellipsoid generalizes; this is another problem addressed by Jacobi.</li>
</ul>
</li>
<li>My interest in this problem was piqued by Jean-Marc Baillard. I put him onto Jacobi's solution without having looked at it in detail myself; and he quickly implemented the solution for an HP-41 calculator(!) which is posted <a href="http://hp41programs.yolasite.com/geod3axial.php">here</a>.</li>
<li>I do not give full citations of the papers here. You can find these in the <a href="http://geographiclib.sourceforge.net/geodesic-papers/biblio.html">Online Geodesic Bibliography</a>; this includes links to online versions of the papers.</li>
<li>An alternative to exploring geodesics using Jacobi's solution is to integrate the equations for the geodesics directly. This is the approach taken by <a href="http://www.math.harvard.edu/~knill/caustic/">Oliver Knill and Michael Teodorescu</a>. However it is difficult to ensure that the long time behavior is correctly modeled with such an approach.</li>
<li>At this point, I have no plans to add the solution of triaxial geodesic problem to <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a>.</li>
<li>If you only want to learn about geodesics on a biaxial ellipsoid (an ellipsoid of revolution), then see <a class="el" href="geodesic.html">Geodesics on an ellipsoid of revolution</a> or the paper<ul>
<li>C. F. F. Karney, <a href="https://dx.doi.org/10.1007/s00190-012-0578-z">Algorithms for geodesics</a>, J. Geodesy 87(1), 43&ndash;55 (2013); DOI: <a href="https://dx.doi.org/10.1007/s00190-012-0578-z">10.1007/s00190-012-0578-z</a>; addenda: <a href="http://geographiclib.sf.net/geod-addenda.html">geod-addenda.html</a>.</li>
</ul>
</li>
</ol>
<h1><a class="anchor" id="triaxial-coords"></a>
Triaxial coordinate systems</h1>
<p>Consider the ellipsoid defined by </p><p class="formulaDsp">
\[ f = \frac{X^2}{a^2} + \frac{Y^2}{b^2} + \frac{Z^2}{c^2} = 1, \]
</p>
<p> 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 <em>geographical</em> latitude and longitude \((\phi, \lambda)\) are defined by </p><p class="formulaDsp">
\[ \frac{\nabla f}{\left| \nabla f\right|} = \left( \begin{array}{c} \cos\phi \cos\lambda \\ \cos\phi \sin\lambda \\ \sin\phi \end{array}\right). \]
</p>
<p> The <em>parametric</em> latitude and longitude \((\phi&#39;, \lambda&#39;)\) are defined by </p><p class="formulaDsp">
\[ \begin{align} X &amp;= a \cos\phi&#39; \cos\lambda&#39;, \\ Y &amp;= b \cos\phi&#39; \sin\lambda&#39;, \\ Z &amp;= c \sin\phi&#39;. \end{align} \]
</p>
<p> Jacobi employed the <em>ellipsoidal</em> latitude and longitude \((\beta, \omega)\) defined by </p><p class="formulaDsp">
\[ \begin{align} X &amp;= a \cos\omega \frac{\sqrt{a^2 - b^2\sin^2\beta - c^2\cos^2\beta}} {\sqrt{a^2 - c^2}}, \\ Y &amp;= b \cos\beta \sin\omega, \\ Z &amp;= c \sin\beta \frac{\sqrt{a^2\sin^2\omega + b^2\cos^2\omega - c^2}} {\sqrt{a^2 - c^2}}. \end{align} \]
</p>
<p> Grid lines of constant \(\beta\) and \(\omega\) are given in Fig. 1.</p>
<center> <div class="image">
<img src="https://upload.wikimedia.org/wikipedia/commons/b/bd/Triaxial_ellipsoid_coordinate_system.svg"  width="419" height="397" alt="Geodesic grid on a triaxial ellipsoid"/>
</div>
<p> Fig. 1 </center><p><br />
Fig. 1: The ellipsoidal grid. The blue (resp. green) lines are lines of constant \(\beta\) (resp. \(\omega\)); the grid spacing is 10&deg;. 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.</p>
<p>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.</p>
<p>The element of length on the ellipsoid in ellipsoidal coordinates is given by </p><p class="formulaDsp">
\[ \begin{align} \frac{ds^2}{(a^2-b^2)\sin^2\omega+(b^2-c^2)\cos^2\beta} &amp;= \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 \\ &amp;\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} \]
</p>
<p>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 </p><p class="formulaDsp">
\[ \begin{align} \omega &amp; \rightarrow -\omega,\\ \beta &amp; \rightarrow \pi-\beta,\\ \alpha &amp; \rightarrow \pi+\alpha, \end{align} \]
</p>
<p> where \(\alpha\) is the heading of a path, relative to a line of constant \(\omega\).</p>
<p>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)\).</p>
<p>Inter-conversions between the three different latitudes and longitudes and the cartesian coordinates are simple algebraic exercises.</p>
<h1><a class="anchor" id="triaxial-jacobi"></a>
Jacobi's solution</h1>
<p>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.</p>
<p>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 <a href="https://books.google.com/books?id=_09tAAAAMAAJ&amp;pg=PA385">letter</a> he wrote to his friend and neighbor Bessel: </p><blockquote class="doxtable">
<p>The day before yesterday, I reduced to quadrature the problem of geodesic lines on an <em>ellipsoid with three unequal axes</em>. They are the simplest formulas in the world, Abelian integrals, which become the well known elliptic integrals if 2 axes are set equal.<br />
K&ouml;nigsberg, 28th Dec. '38. </p>
</blockquote>
<p>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).</p>
<p>Here is the solution, exactly as given by Jacobi <a href="https://books.google.com/books?id=Rh8GAAAAYAAJ&amp;pg=PA268">here</a> (with minor changes in notation): </p><p class="formulaDsp">
\[ \begin{align} \delta &amp;= \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}}\\ &amp;\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} \]
</p>
<p> As Jacobi notes &quot;a function of the angle \(\beta\) equals a function of the angle \(\omega\). These two functions are just Abelian integrals&hellip;&quot; 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 </p><p class="formulaDsp">
\[ \gamma = (b^2-c^2)\cos^2\beta\sin^2\alpha-(a^2-b^2)\sin^2\omega\cos^2\alpha \]
</p>
<p> 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) <a href="https://books.google.com/books?id=hGMSAAAAIAAJ&amp;pg=PA9">&sect;&sect;583&ndash;584</a> 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 </p><p class="formulaDsp">
\[ \begin{align} \frac{ds}{(b^2-c^2)\cos^2\beta + (a^2-b^2)\sin^2\omega} &amp;= \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}}\\ &amp;= \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} \]
</p>
<p> An alternative expression for the distance is </p><p class="formulaDsp">
\[ \begin{align} ds &amp;= \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}}\\ &amp;\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} \]
</p>
<p>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 integrals are periodic, so the behavior of an arbitrarily long geodesic is entirely captured by tabulating the integrals over a single period.</p>
<p>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.</p>
<h1><a class="anchor" id="triaxial-survey"></a>
Survey of triaxial geodesics</h1>
<p>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.</p>
<p>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&ndash;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.</p>
<p>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.</p>
<center> <div class="image">
<img src="https://upload.wikimedia.org/wikipedia/commons/a/a9/Circumpolar_geodesic_on_a_triaxial_ellipsoid_case_A.svg"  width="419" height="397" alt="Example of a circumpolar geodesic on a
triaxial ellipsoid"/>
</div>
<p> Fig. 2 </center><p><br />
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\).</p>
<center> <div class="image">
<img src="https://upload.wikimedia.org/wikipedia/commons/f/f3/Circumpolar_geodesic_on_a_triaxial_ellipsoid_case_B.svg"  width="419" height="397" alt="Another example of a circumpolar geodesic on a
triaxial ellipsoid"/>
</div>
<p> Fig. 3 </center><p><br />
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\).</p>
<p>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.</p>
<center> <div class="image">
<img src="https://upload.wikimedia.org/wikipedia/commons/8/83/Transpolar_geodesic_on_a_triaxial_ellipsoid_case_A.svg"  width="419" height="397" alt="Example of a transpolar geodesic on a
triaxial ellipsoid"/>
</div>
<p> Fig. 4 </center><p><br />
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\).</p>
<center> <div class="image">
<img src="https://upload.wikimedia.org/wikipedia/commons/9/9c/Transpolar_geodesic_on_a_triaxial_ellipsoid_case_B.svg"  width="419" height="397" alt="Another example of a transpolar geodesic on a
triaxial ellipsoid"/>
</div>
<p> Fig. 5 </center><p><br />
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\).</p>
<p>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.</p>
<center> <div class="image">
<img src="https://upload.wikimedia.org/wikipedia/commons/e/e5/Unstable_umbilical_geodesic_on_a_triaxial_ellipsoid.svg"  width="419" height="397" alt="Example of an umbilical geodesic on a
triaxial ellipsoid"/>
</div>
<p> Fig. 6 </center><p><br />
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.</p>
<p>Umbilical geodesics enjoy several interesting properties.</p><ul>
<li>Through any point on the ellipsoid, there are two umbilical geodesics.</li>
<li>The geodesic distance between opposite umbilical points is the same regardless of the initial direction of the geodesic.</li>
<li>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 <em>unstable</em>. 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.).</li>
</ul>
<h1><a class="anchor" id="triaxial-stab"></a>
The stability of closed geodesics</h1>
<p>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 </p><p class="formulaDsp">
\[ \frac {d^2m}{ds^2} + Km = 0 \]
</p>
<p> 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 <a href="http://dlmf.nist.gov/28.29">DLMF, &sect;28.29</a>. Using the notation of &sect;3 of <a href="https://dx.doi.org/10.1007/s00190-012-0578-z">Algorithms for geodesics</a>, 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 </p><p class="formulaDsp">
\[ {\cal M} = \left(\begin{array}{cc} M_{12} &amp; m_{12}\\ -\frac{1 - M_{12}M_{21}}{m_{12}} &amp; M_{21} \end{array}\right). \]
</p>
<p> Because \(\mathrm{det}\,{\cal M} = 1\), the eigenvalues are determined by \(\mathrm{tr}\,{\cal M}\). In particular if \(\left|\mathrm{tr}\,{\cal M}\right| &lt; 2\), we have \(\left|\lambda_{1,2}\right| = 1\) and the solution is stable; if \(\left|\mathrm{tr}\,{\cal M}\right| &gt; 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.</p>
<p>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.)</p>
<h1><a class="anchor" id="triaxial-inverse"></a>
The inverse problem</h1>
<p>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)\).</p><ul>
<li>If the point is an umbilical point, all the lines meet at the opposite umbilical point.</li>
<li>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 <a href="https://dx.doi.org/10.1080/10586458.2003.10504515">Sinclair (2003)</a>.</li>
<li>All geodesics intersect (or, in the case of \(\alpha_1 = 0\) or \(\pi\), touch) the line \(\omega = \omega_1 + \pi\).</li>
<li>All geodesics intersect (or, in the case of \(\alpha_1 = \pm\pi/2\), touch) the line \(\beta = -\beta_1\).</li>
<li>Two geodesics with azimuths \(\pm\alpha_1\) first intersect on \(\omega = \omega_1 + \pi\) and their lengths to the point of intersection are equal.</li>
<li>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.</li>
</ul>
<p>(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\).</p><ul>
<li>The geodesics completely fill the portion of the ellipsoid satisfying \(\beta \le -\beta_1\).</li>
<li>None of geodesics intersect any other geodesics.</li>
<li>Any initial portion of these geodesics is a shortest path.</li>
<li>Each geodesic intersects the line \(\beta = \beta_2\), where \(\beta_1 &lt; \beta_2 &lt; -\beta_1\), exactly once.</li>
<li>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\).</li>
<li>If \(\beta_2 = \pm \beta_1\), then the previous two assertions need to be modified similarly to the case for an ellipsoid of revolution.</li>
</ul>
<p>These properties show that the inverse problem can be solved using techniques similar to those employed for an ellipsoid of revolution (see &sect;4 of <a href="https://dx.doi.org/10.1007/s00190-012-0578-z">Algorithms for geodesics</a>).</p><ul>
<li>If the points are opposite umbilical points, an arbitrary \(\alpha_1\) may be chosen.</li>
<li>If the points are neighboring umbilical points, the shortest path lies on the ellipse \(y = 0\).</li>
<li>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\).</li>
<li>Treat the cases where the geodesic might follow a line of constant \(\beta\). There are two such cases: (a)&#160;the points lie on the ellipse \(z = 0\) on a general ellipsoid and (b)&#160;the points lie on an ellipse whose major axis is the \(x\) axis on a prolate ellipsoid ( \(a = b &gt; 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).</li>
<li>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 <em>hybrid</em> 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}\).</li>
</ul>
<p>The shortest path found by this method is unique unless:</p><ul>
<li>The length of the geodesic vanishes \(s_{12}=0\), in which case any constant can be added to the azimuths.</li>
<li>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 <em>forward</em> azimuth at the second umbilical point, i.e., its azimuth immediately <em>after</em> passage through the umbilical point.</li>
<li>\(\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\).</li>
</ul>
<h1><a class="anchor" id="triaxial-conformal"></a>
Jacobi's conformal projection</h1>
<p>This material is now on its own page; see <a class="el" href="jacobi.html">Jacobi's conformal projection</a>.</p>
<center> Back to <a class="el" href="geodesic.html">Geodesics on an ellipsoid of revolution</a>. Forward to <a class="el" href="jacobi.html">Jacobi's conformal projection</a>. Up to <a class="el" href="index.html#contents">Contents</a>. </center> </div></div><!-- contents -->
<!-- start footer part -->
<hr class="footer"/><address class="footer"><small>
Generated by &#160;<a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/>
</a> 1.8.9.1
</small></address>
</body>
</html>

APT Browse - Built by Thomas Orozco.