/usr/share/doc/geographiclib/html/rhumb.html is in geographiclib-doc 1.49-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 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 | <!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.13"/>
<meta name="viewport" content="width=device-width, initial-scale=1"/>
<title>GeographicLib: Rhumb lines</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 type="text/javascript" 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 id="projectalign" style="padding-left: 0.5em;">
<div id="projectname">GeographicLib
 <span id="projectnumber">1.49</span>
</div>
</td>
</tr>
</tbody>
</table>
</div>
<!-- end header part -->
<!-- Generated by Doxygen 1.8.13 -->
<script type="text/javascript" src="menudata.js"></script>
<script type="text/javascript" src="menu.js"></script>
<script type="text/javascript">
$(function() {
initMenu('',false,false,'search.php','Search');
});
</script>
<div id="main-nav"></div>
</div><!-- top -->
<div class="header">
<div class="headertitle">
<div class="title"><a class="el" href="classGeographicLib_1_1Rhumb.html" title="Solve of the direct and inverse rhumb problems. ">Rhumb</a> lines </div> </div>
</div><!--header-->
<div class="contents">
<div class="textblock"><center> Back to <a class="el" href="jacobi.html">Jacobi's conformal projection</a>. Forward to <a class="el" href="greatellipse.html">Great Ellipses</a>. Up to <a class="el" href="index.html#contents">Contents</a>. </center><p>The <a class="el" href="classGeographicLib_1_1Rhumb.html" title="Solve of the direct and inverse rhumb problems. ">Rhumb</a> and <a class="el" href="classGeographicLib_1_1RhumbLine.html" title="Find a sequence of points on a single rhumb line. ">RhumbLine</a> classes together with the <a href="RhumbSolve.1.html">RhumbSolve</a> utility perform rhumb line calculations. A rhumb line (also called a loxodrome) is a line of constant azimuth on the surface of the ellipsoid. It is important for historical reasons because sailing with a constant compass heading is simpler than sailing the shorter geodesic course; see Osborne (2013). The formulation of the problem on an ellipsoid is covered by Smart (1946) and Carlton-Wippern (1992). Computational approaches are given by Williams (1950) and Bennett (1996). My interest in this problem was piqued by Botnev and Ustinov (2014) who discuss various techniques to improve the accuracy of rhumb line calculations. The items of interest here are:</p><ul>
<li>Review of accurate formulas for the auxiliary latitudes.</li>
<li>The calculation of the area under a rhumb line.</li>
<li>Using divided differences to compute \(\mu_{12}/\psi_{12}\) maintaining full accuracy. Two cases are treated:<ul>
<li>If \(f\) is small, Krüger's series for the transverse Mercator projection relate \(\chi\) and \(\mu\).</li>
<li>For arbitrary \(f\), the divided difference formula for incomplete elliptic integrals of the second relates \(\beta\) and \(\mu\).</li>
</ul>
</li>
<li>Extending Clenshaw summation to compute the divided difference of a trigonometric sum.</li>
</ul>
<p>Go to</p><ul>
<li><a class="el" href="rhumb.html#rhumbform">Formulation of the rhumb line problem</a></li>
<li><a class="el" href="rhumb.html#rhumblat">Determining the auxiliary latitudes</a></li>
<li><a class="el" href="rhumb.html#rhumbarea">The area under a rhumb line</a></li>
<li><a class="el" href="rhumb.html#divideddiffs">Use of divided differences</a></li>
<li><a class="el" href="rhumb.html#dividedclenshaw">Clenshaw evaluation of differenced sums</a></li>
</ul>
<p>References:</p><ul>
<li>G. G. Bennett, <a href="https://doi.org/10.1017/S0373463300013151">Practical Rhumb Line Calculations on the Spheroid</a> J. Navigation 49(1), 112–119 (1996).</li>
<li>F. W. Bessel, <a href="https://doi.org/10.1002/asna.201011352">The calculation of longitude and latitude from geodesic measurements (1825)</a>, Astron. Nachr. 331(8), 852–861 (2010); translated by C. F. F. Karney and R. E. Deakin; preprint: <a href="https://arxiv.org/abs/0908.1824">arXiv:0908.1824</a>.</li>
<li>V.A. Botnev, S.M. Ustinov, <a href="http://ntv.spbstu.ru/fulltext/T3.198.2014_05.PDF">Metody resheniya pryamoy i obratnoy geodezicheskikh zadach s vysokoy tochnost'yu </a> (Methods for direct and inverse geodesic problems solving with high precision), St. Petersburg State Polytechnical University Journal 3(198), 49–58 (2014).</li>
<li>K. C. Carlton-Wippern <a href="https://doi.org/10.1017/S0373463300010791">On Loxodromic Navigation</a>, J. Navigation 45(2), 292–297 (1992).</li>
<li>K. E. Engsager and K. Poder, <a href="http://icaci.org/files/documents/ICC_proceedings/ICC2007/documents/doc/THEME 2/oral 1/2.1.2 A HIGHLY ACCURATE WORLD WIDE ALGORITHM FOR THE TRANSVE.doc">A highly accurate world wide algorithm for the transverse Mercator mapping (almost)</a>, Proc. XXIII Intl. Cartographic Conf. (ICC2007), Moscow (2007).</li>
<li>F. R. Helmert, <a href="https://doi.org/10.5281/zenodo.32050">Mathematical and Physical Theories of Higher Geodesy, Part 1 (1880)</a>, Aeronautical Chart and Information Center (St. Louis, 1964), Chaps. 5–7.</li>
<li>W. M. Kahan and R. J. Fateman, <a href="https://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf">Symbolic computation of divided differences</a>, SIGSAM Bull. 33(3), 7–28 (1999) DOI: <a href="https://doi.org/10.1145/334714.334716">10.1145/334714.334716</a>.</li>
<li>C. F. F. Karney, <a href="https://doi.org/10.1007/s00190-011-0445-3">Transverse Mercator with an accuracy of a few nanometers</a>, J. Geodesy 85(8), 475–485 (Aug. 2011); addenda: <a href="https://geographiclib.sourceforge.io/tm-addenda.html">tm-addenda.html</a>; preprint: <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>; resource page: <a href="https://geographiclib.sourceforge.io/tm.html">tm.html</a>.</li>
<li>C. F. F. Karney, <a href="https://doi.org/10.1007/s00190-012-0578-z">Algorithms for geodesics</a>, J. Geodesy 87(1), 43–55 (2013); DOI: <a href="https://doi.org/10.1007/s00190-012-0578-z">10.1007/s00190-012-0578-z</a>; addenda: <a href="https://geographiclib.sourceforge.io/geod-addenda.html">geod-addenda.html</a>; resource page: <a href="https://geographiclib.sourceforge.io/geod.html">geod.html</a>.</li>
<li>L. Krüger, <a href="https://doi.org/10.2312/GFZ.b103-krueger28">Konforme Abbildung des Erdellipsoids in der Ebene</a> (Conformal mapping of the ellipsoidal earth to the plane), Royal Prussian Geodetic Institute, New Series 52, 172 pp. (1912).</li>
<li>P. Osborne, <a href="https://doi.org/10.5281/zenodo.35392">The Mercator Projections</a> (2013), §2.5 and §6.5; DOI: <a href="https://doi.org/10.5281/zenodo.35392">10.5281/zenodo.35392</a>.</li>
<li>W. M. Smart <a href="https://doi.org/10.1093/mnras/106.2.124">On a Problem in Navigation</a> MNRAS 106(2), 124–127 (1946).</li>
<li>J. E. D. Williams <a href="https://doi.org/10.1017/S0373463300045549">Loxodromic Distances on the Terrestrial Spheroid Journal</a>, J. Navigation 3(2), 133–140 (1950)</li>
</ul>
<h1><a class="anchor" id="rhumbform"></a>
Formulation of the rhumb line problem</h1>
<p>The rhumb line can formulated in terms of three types of latitude, the isometric latitude \(\psi\) (this serves to define the course of rhumb lines), the rectifying latitude \(\mu\) (needed for computing distances along rhumb lines), and the parametric latitude \(\beta\) (needed for dealing with rhumb lines that run along a parallel). These are defined in terms of the geographical latitude \(\phi\) by </p><p class="formulaDsp">
\[ \begin{align} \frac{d\psi}{d\phi} &= \frac{\rho}R, \\ \frac{d\mu}{d\phi} &= \frac{\pi}{2M}\rho, \\ \end{align} \]
</p>
<p> where \(\rho\) is the meridional radius of curvature, \(R = a\cos\beta\) is the radius of a circle of latitude, and \(M\) is the length of a quarter meridian (from the equator to the pole).</p>
<p><a class="el" href="classGeographicLib_1_1Rhumb.html" title="Solve of the direct and inverse rhumb problems. ">Rhumb</a> lines are straight in the Mercator projection, which maps a point with geographical coordinates \( (\phi,\lambda) \) on the ellipsoid to a point \( (\psi,\lambda) \). The azimuth \(\alpha_{12}\) of a rhumb line from \( (\phi_1,\lambda_1) \) to \( (\phi_2,\lambda_2) \) is thus given by </p><p class="formulaDsp">
\[ \tan\alpha_{12} = \frac{\lambda_{12}}{\psi_{12}}, \]
</p>
<p> where the quadrant of \(\alpha_{12}\) is determined by the signs of the numerator and denominator of the right-hand side, \(\lambda_{12} = \lambda_2 - \lambda_1\), \(\psi_{12} = \psi_2 - \psi_1\), and, typically, \(\lambda_{12}\) is reduced to the range \([-\pi,\pi]\) (thus giving the course of the <em>shortest</em> rhumb line).</p>
<p>The distance is given by </p><p class="formulaDsp">
\[ \begin{align} s_{12} &= \frac {2M}{\pi} \mu_{12} \sec\alpha_{12} \\ &= \frac {2M}{\pi} \frac{\mu_{12}}{\psi_{12}} \sqrt{\lambda_{12}^2 + \psi_{12}^2}, \end{align} \]
</p>
<p> where \(\mu_{12} = \mu_2 - \mu_1\). This relation is indeterminate if \(\phi_1 = \phi_2\), so we take the limits \(\psi_{12}\rightarrow0\) and \(\mu_{12}\rightarrow0\) and apply L'Hôpital's rule to yield </p><p class="formulaDsp">
\[ \begin{align} s_{12} &= \frac {2M}{\pi} \frac{d\mu}{d\psi} \left|\lambda_{12}\right|\\ &=a \cos\beta \left|\lambda_{12}\right|, \end{align} \]
</p>
<p> where the last relation is given by the defining equations for \(\psi\) and \(\mu\).</p>
<p>This provides a complete solution for rhumb lines. This formulation entirely encapsulates the ellipsoidal shape of the earth via the auxiliary latitudes; thus in the <a class="el" href="classGeographicLib_1_1Rhumb.html" title="Solve of the direct and inverse rhumb problems. ">Rhumb</a> class, the ellipsoidal generalization is handled by the <a class="el" href="classGeographicLib_1_1Ellipsoid.html" title="Properties of an ellipsoid. ">Ellipsoid</a> class where the auxiliary latitudes are defined.</p>
<h1><a class="anchor" id="rhumblat"></a>
Determining the auxiliary latitudes</h1>
<p>Here we brief develop the necessary formulas for the auxiliary latitudes.</p>
<p><b>Isometric latitude:</b> The equation for \(\psi\) can be integrated to give </p><p class="formulaDsp">
\[ \psi = \sinh^{-1}\tan\phi - e\tanh^{-1}(e \sin\phi), \]
</p>
<p> where \(e = \sqrt{f(2-f)}\) is the eccentricity of the ellipsoid. To invert this equation (to give \(\phi\) in terms of \(\psi\)), it is convenient to introduce the variables \(\tau=\tan\phi\) and \(\tau' = \tan\chi = \sinh\psi\) ( \(\chi\) is the conformal latitude) which are related by (Karney, 2011) </p><p class="formulaDsp">
\[ \begin{align} \tau' &= \tau \sqrt{1 + \sigma^2} - \sigma \sqrt{1 + \tau^2}, \\ \sigma &= \sinh\bigl(e \tanh^{-1}(e \tau/\sqrt{1 + \tau^2}) \bigr). \end{align} \]
</p>
<p> The equation for \(\tau'\) can be inverted to give \(\tau\) in terms of \(\tau'\) using Newton's method with \(\tau_0 = \tau'/(1-e^2)\) as a starting guess; and, of course, \(\phi\) and \(\psi\) are then given by </p><p class="formulaDsp">
\[ \begin{align} \phi &= \tan^{-1}\tau,\\ \psi &= \sinh^{-1}\tau'. \end{align} \]
</p>
<p> This allows conversions to and from \(\psi\) to be carried out for any value of the flattening \(f\); these conversions are implemented by <a class="el" href="classGeographicLib_1_1Ellipsoid.html#a35c553ad95bb70e6ecf154f3cfb5e7b8">Ellipsoid::IsometricLatitude</a> and <a class="el" href="classGeographicLib_1_1Ellipsoid.html#a6e549871bae33dd72b6947077eec9277">Ellipsoid::InverseIsometricLatitude</a>. (For prolate ellipsoids, \(f\) is negative and \(e\) is imaginary, and the equation for \(\sigma\) needs to be recast in terms of the tangent function.)</p>
<p>For small values of \(f\), Engsager and Poder (2007) express \(\chi\) as a trigonometric series in \(\phi\). This series can be reverted to give a series for \(\phi\) in terms of \(\chi\). Both series can be efficiently evaluated using Clenshaw summation and this provides a fast non-iterative way of making the conversions.</p>
<p><b>Parametric latitude:</b> This is given by </p><p class="formulaDsp">
\[ \tan\beta = (1-f)\tan\phi, \]
</p>
<p> which allows rapid and accurate conversions; these conversions are implemented by <a class="el" href="classGeographicLib_1_1Ellipsoid.html#ad8041795c821fa6247db7ef6e5bb5c6a">Ellipsoid::ParametricLatitude</a> and <a class="el" href="classGeographicLib_1_1Ellipsoid.html#a7218afae76946bf4094c65a9382e2263">Ellipsoid::InverseParametricLatitude</a>.</p>
<p><b>Rectifying latitude:</b> Solving for distances on the meridian is naturally carried out in terms of the parametric latitude instead of the geographical latitude. This leads to a simpler elliptic integral which is easier to evaluate and to invert. Helmert (1880, §5.11) notes that the series for meridian distance in terms of \(\beta\) converge faster than the corresponding ones in terms of \(\phi\).</p>
<p>In terms of \(\beta\), the rectifying latitude is given by </p><p class="formulaDsp">
\[ \mu = \frac{\pi}{2E(ie')} E(\beta,ie'), \]
</p>
<p> where \(e'=e/\sqrt{1-e^2}\) is the second eccentricity and \(E(k)\) and \(E(\phi,k)\) are the complete and incomplete elliptic integrals of the second kind; see <a href="http://dlmf.nist.gov/19.2.ii">http://dlmf.nist.gov/19.2.ii</a>. These can be evaluated accurately for arbitrary flattening using the method given in <a href="http://dlmf.nist.gov/19.36.i">http://dlmf.nist.gov/19.36.i</a>. To find \(\beta\) in terms of \(\mu\), Newton's method can be used (noting that \(dE(\phi,k)/d\phi = \sqrt{1 - k^2\sin^2\phi}\)); for a starting guess use \(\beta_0 = \mu\) (or use the first terms in the reverted series; see below). These conversions are implemented by <a class="el" href="classGeographicLib_1_1Ellipsoid.html#a31fe44305d950a8ceb528d33de15f183">Ellipsoid::RectifyingLatitude</a> and <a class="el" href="classGeographicLib_1_1Ellipsoid.html#a76442fb6a136f514f61386c3f237777b">Ellipsoid::InverseRectifyingLatitude</a>.</p>
<p>If the flattening is small, \(\mu\) can be expressed as a series in various ways. The most economical series is in terms of the third flattening \( n = f/(2-f)\) and was found by Bessel (1825); see Eq. 5.5.7 of Helmert (1880). Helmert (1880, Eq. 5.6.8) also gives the reverted series (finding \(\beta\) given \(\mu\)). These series are used by the <a class="el" href="classGeographicLib_1_1Geodesic.html" title="Geodesic calculations ">Geodesic</a> class where the coefficients are \(C_{1j}\) and \(C_{1j}'\); see Eqs. (18) and (21) of Karney (2013) and <a class="el" href="geodesic.html#geodseries">Expansions for geodesics</a>. (Make the replacements, \(\sigma\rightarrow\beta\), \(\tau\rightarrow\mu\), \(\epsilon\rightarrow n\), to convert the notation of the geodesic problem to the problem at hand.) The series are evaluated by Clenshaw summation.</p>
<p>These relations allow inter-conversions between the various latitudes \(\phi\), \(\psi\), \(\chi\), \(\beta\), \(\mu\) to be carried out simply and accurately. The approaches using series apply only if \(f\) is small. The others apply for arbitrary values of \(f\).</p>
<h1><a class="anchor" id="rhumbarea"></a>
The area under a rhumb line</h1>
<p>The area between a rhumb line and the equator is given by </p><p class="formulaDsp">
\[ S_{12} = \int_{\lambda_1}^{\lambda_2} c^2 \sin\xi \,d\lambda, \]
</p>
<p> where \(c\) is the authalic radius and \(\xi\) is the authalic latitude. Express \(\sin\xi\) in terms of the conformal latitude \(\chi\) and expand as a series in the third flattening \(n\). This can be expressed in terms of powers of \(\sin\chi\). Substitute \(\sin\chi=\tanh\psi\), where \(\psi\) is the isometric latitude. For a rhumb line we have </p><p class="formulaDsp">
\[ \begin{align} \psi &= m(\lambda-\lambda_0),\\ m &= \frac{\psi_2 - \psi_1}{\lambda_2 - \lambda_1}. \end{align} \]
</p>
<p> Performing the integral over \(\lambda\) gives </p><p class="formulaDsp">
\[ S_{12} = c^2 \lambda_{12} \left<\sin\xi\right>_{12}, \]
</p>
<p> where \(\left<\sin\xi\right>_{12}\) is the mean value of \(\sin\xi\) given by </p><p class="formulaDsp">
\[ \begin{align} \left<\sin\xi\right>_{12} &= \frac{S(\chi_2) - S(\chi_1)}{\psi_2 - \psi_1},\\ S(\chi) &= \log\sec\chi + \sum_{l=1} R_l\cos(2l\chi), \end{align} \]
</p>
<p> \(\log\) is the natural logarithm, and \(R_l = O(n^l)\) is given as a series in \(n\) below. In the spherical limit, the sum vanishes.</p>
<p>Note the simple way that longitude enters into the expression for the area. In the limit \(\chi_2 \rightarrow \chi_1\), we can apply l'Hôpital's rule </p><p class="formulaDsp">
\[ \left<\sin\xi\right>_{12} \rightarrow \frac{dS(\chi_1)/d\chi_1}{\sec\chi_1} =\sin\xi_1, \]
</p>
<p> as expected.</p>
<p>In order to maintain accuracy, particularly for rhumb lines which nearly follow a parallel, evaluate \(\left<\sin\xi\right>_{12}\) using divided differences (see the <a class="el" href="rhumb.html#divideddiffs">next section</a>) and use Clenshaw summation to evaluate the sum (see the <a class="el" href="rhumb.html#dividedclenshaw">last section</a>).</p>
<p>Here is the series expansion accurate to 10th order, found by the Maxima script <a href="rhumbarea.mac">rhumbarea.mac</a>:</p>
<pre class="fragment">R[1] = - 1/3 * n
+ 22/45 * n^2
- 356/945 * n^3
+ 1772/14175 * n^4
+ 41662/467775 * n^5
- 114456994/638512875 * n^6
+ 258618446/1915538625 * n^7
- 1053168268/37574026875 * n^8
- 9127715873002/194896477400625 * n^9
+ 33380126058386/656284056553125 * n^10;
R[2] = - 2/15 * n^2
+ 106/315 * n^3
- 1747/4725 * n^4
+ 18118/155925 * n^5
+ 51304574/212837625 * n^6
- 248174686/638512875 * n^7
+ 2800191349/14801889375 * n^8
+ 10890707749202/64965492466875 * n^9
- 3594078400868794/10719306257034375 * n^10;
R[3] = - 31/315 * n^3
+ 104/315 * n^4
- 23011/51975 * n^5
+ 1554472/14189175 * n^6
+ 114450437/212837625 * n^7
- 8934064508/10854718875 * n^8
+ 4913033737121/21655164155625 * n^9
+ 591251098891888/714620417135625 * n^10;
R[4] = - 41/420 * n^4
+ 274/693 * n^5
- 1228489/2027025 * n^6
+ 3861434/42567525 * n^7
+ 1788295991/1550674125 * n^8
- 215233237178/123743795175 * n^9
+ 95577582133463/714620417135625 * n^10;
R[5] = - 668/5775 * n^5
+ 1092376/2027025 * n^6
- 3966679/4343625 * n^7
+ 359094172/10854718875 * n^8
+ 7597613999411/3093594879375 * n^9
- 378396252233936/102088631019375 * n^10;
R[6] = - 313076/2027025 * n^6
+ 4892722/6081075 * n^7
- 1234918799/834978375 * n^8
- 74958999806/618718975875 * n^9
+ 48696857431916/9280784638125 * n^10;
R[7] = - 3189007/14189175 * n^7
+ 930092876/723647925 * n^8
- 522477774212/206239658625 * n^9
- 2163049830386/4331032831125 * n^10;
R[8] = - 673429061/1929727800 * n^8
+ 16523158892/7638505875 * n^9
- 85076917909/18749059875 * n^10;
R[9] = - 39191022457/68746552875 * n^9
+ 260863656866/68746552875 * n^10;
R[10] = - 22228737368/22915517625 * n^10;
</pre><h1><a class="anchor" id="divideddiffs"></a>
Use of divided differences</h1>
<p>Despite our ability to compute latitudes accurately, the way that distances enter into the solution involves the ratio \(\mu_{12}/\psi_{12}\); the numerical calculation of this term is subject to catastrophic round-off errors when \(\phi_1\) and \(\phi_2\) are close. A simple solution, suggested by Bennett (1996), is to extend the treatment of rhumb lines along parallels to this case, i.e., to replace the ratio by the derivative evaluated at the midpoint. However this is not entirely satisfactory: you have to pick the transition point where the derivative takes over from the ratio of differences; and, near this transition point, many bits of accuracy will be lost (I estimate that about 1/3 of bits are lost, leading to errors on the order of 0.1 mm for double precision). Note too that this problem crops up even in the spherical limit.</p>
<p>It turns out that we can do substantially better than this and maintain full double precision accuracy. Indeed Botnev and Ustinov provide formulas which allow \(\mu_{12}/\psi_{12}\) to be computed accurately (but the formula for \(\mu_{12}\) applies only for small flattening).</p>
<p>Here I give a more systematic treatment using the algebra of <em>divided differences</em>. Many readers will be already using this technique, e.g., when writing, for example, </p><p class="formulaDsp">
\[ \frac{\sin(2x) - \sin(2y)}{x-y} = 2\frac{\sin(x-y)}{x-y}\cos(x+y). \]
</p>
<p> However, Kahan and Fateman (1999) provide an extensive set of rules which allow the technique to be applied to a wide range of problems. (The classes <a class="el" href="classGeographicLib_1_1LambertConformalConic.html" title="Lambert conformal conic projection. ">LambertConformalConic</a> and <a class="el" href="classGeographicLib_1_1AlbersEqualArea.html" title="Albers equal area conic projection. ">AlbersEqualArea</a> use this technique to improve the accuracy.)</p>
<p>To illustrate the technique, consider the relation between \(\psi\) and \(\chi\), </p><p class="formulaDsp">
\[ \psi = \sinh^{-1}\tan\chi. \]
</p>
<p> The divided difference operator is defined by </p><p class="formulaDsp">
\[ \Delta[f](x,y) = \begin{cases} df(x)/dx, & \text{if $x=y$,} \\ \displaystyle\frac{f(x)-f(y)}{x-y}, & \text{otherwise.} \end{cases} \]
</p>
<p> Many of the rules for differentiation apply to divided differences; in particular, we can use the chain rule to write </p><p class="formulaDsp">
\[ \frac{\psi_1 - \psi_2}{\chi_1 - \chi_2} = \Delta[\sinh^{-1}](\tan\chi_1,\tan\chi_2) \times \Delta[\tan](\chi_1,\chi_2). \]
</p>
<p>Kahan and Fateman catalog the divided difference formulas for all the elementary functions. In this case, we need </p><p class="formulaDsp">
\[ \begin{align} \Delta[\tan](x,y) &= \frac{\tan(x-y)}{x-y} (1 + \tan x\tan y),\\ \Delta[\sinh^{-1}](x,y) &= \frac1{x-y} \sinh^{-1}\biggl(\frac{(x-y)(x+y)}{x\sqrt{1+y^2}+y\sqrt{1+x^2}}\biggr). \end{align} \]
</p>
<p> The crucial point in these formulas is that the right hand sides can be evaluated accurately even if \(x-y\) is small. (I've only included the basic results here; Kahan and Fateman supplement these with the derivatives if \(x-y\) vanishes or the direct ratios if \(x-y\) is not small.)</p>
<p>To complete the computation of \(\mu_{12}/\psi_{12}\), we now need \((\mu_1 - \mu_2)/(\chi_1 - \chi_2)\), i.e., we need the divided difference of the function that converts the conformal latitude to rectifying latitude. We could go through the chain of relations between these quantities. However, we can take a short cut and recognize that this function was given as a trigonometric series by Krüger (1912) in his development of the transverse Mercator projection. This is also given by Eqs. (5) and (35) of Karney (2011) with the replacements \(\zeta \rightarrow \mu\) and \(\zeta' \rightarrow \chi\). The coefficients appearing in this series and the reverted series Eqs. (6) and (36) are already computed by the <a class="el" href="classGeographicLib_1_1TransverseMercator.html" title="Transverse Mercator projection. ">TransverseMercator</a> class. The series can be readily converted to divided difference form (see the example at the beginning of this section) and Clenshaw summation can be used to evaluate it (see below).</p>
<p>The approach of using the series for the transverse Mercator projection limits the applicability of the method to \(\left|f\right|\lt 0.01\). If we want to extend the method to arbitrary flattening we need to compute \(\Delta[E](x,y;k)\). The necessary relation is the "addition
theorem" for the incomplete elliptic integral of the second kind given in <a href="http://dlmf.nist.gov/19.11.E2">http://dlmf.nist.gov/19.11.E2</a>. This can be converted in the following divided difference formula </p><p class="formulaDsp">
\[ \Delta[E](x,y;k) =\begin{cases} \sqrt{1 - k^2\sin^2x}, & \text{if $x=y$,} \\ \displaystyle \frac{E(x,k)-E(y,k)}{x-y}, & \text{if $xy \le0$,}\\ \displaystyle \biggl(\frac{E(z,k)}{\sin z} - k^2 \sin x \sin y\biggr)\frac{\sin z}{x-y}, &\text{otherwise,} \end{cases} \]
</p>
<p> where the angle \(z\) is given by </p><p class="formulaDsp">
\[ \begin{align} \sin z &= \frac{2t}{1 + t^2},\quad\cos z = \frac{(1-t)(1+t)}{1 + t^2},\\ t &= \frac{(x-y)\Delta[\sin](x,y)} {\sin x\sqrt{1 - k^2\sin^2y} + \sin y\sqrt{1 - k^2\sin^2x}} \frac{\sin x + \sin y}{\cos x + \cos y}. \end{align} \]
</p>
<p> We also need to apply the divided difference formulas to the conversions from \(\phi\) to \(\beta\) and \(\phi\) to \(\psi\); but these all involve elementary functions and the techniques given in Kahan and Fateman can be used.</p>
<p>The end result is that the <a class="el" href="classGeographicLib_1_1Rhumb.html" title="Solve of the direct and inverse rhumb problems. ">Rhumb</a> class allows the computation of all rhumb lines for any flattening with full double precision accuracy (the maximum error is about 10 nanometers). I've kept the implementation simple, which results in rather a lot of repeated evaluations of quantities. However, in this case, the need for clarity trumps the desire for speed.</p>
<h1><a class="anchor" id="dividedclenshaw"></a>
Clenshaw evaluation of differenced sums</h1>
<p>The use of <a href="https://en.wikipedia.org/wiki/Clenshaw_algorithm#Meridian_arc_length_on_the_ellipsoid">Clenshaw summation</a> for summing series of the form, </p><p class="formulaDsp">
\[ g(x) = \sum_{k=1}^N c_k \sin kx, \]
</p>
<p> is well established. However when computing divided differences, we are interested in evaluating </p><p class="formulaDsp">
\[ g(x)-g(y) = \sum_{k=1}^N 2 c_k \sin\bigl({\textstyle\frac12}k(x-y)\bigr) \cos\bigl({\textstyle\frac12}k(x+y)\bigr). \]
</p>
<p> Clenshaw summation can be used in this case if we simultaneously compute \(g(x)+g(y)\) and perform a matrix summation, </p><p class="formulaDsp">
\[ \mathsf G(x,y) = \begin{bmatrix} (g(x) + g(y)) / 2\\ (g(x) - g(y)) / (x - y) \end{bmatrix} = \sum_{k=1}^N c_k \mathsf F_k(x,y), \]
</p>
<p> where </p><p class="formulaDsp">
\[ \mathsf F_k(x,y)= \begin{bmatrix} \cos kd \sin kp\\ {\displaystyle\frac{\sin kd}d} \cos kp \end{bmatrix}, \]
</p>
<p> \(d=\frac12(x-y)\), \(p=\frac12(x+y)\), and, in the limit \(d\rightarrow0\), \((\sin kd)/d \rightarrow k\). The first element of \(\mathsf G(x,y)\) is the average value of \(g\) and the second element, the divided difference, is the average slope. \(\mathsf F_k(x,y)\) satisfies the recurrence relation </p><p class="formulaDsp">
\[ \mathsf F_{k+1}(x,y) = \mathsf A(x,y) \mathsf F_k(x,y) - \mathsf F_{k-1}(x,y), \]
</p>
<p> where </p><p class="formulaDsp">
\[ \mathsf A(x,y) = 2\begin{bmatrix} \cos d \cos p & -d\sin d \sin p \\ - {\displaystyle\frac{\sin d}d} \sin p & \cos d \cos p \end{bmatrix}, \]
</p>
<p> and \(\lim_{d\rightarrow0}(\sin d)/d = 1\). The standard Clenshaw algorithm can now be applied to yield </p><p class="formulaDsp">
\[ \begin{align} \mathsf B_{N+1} &= \mathsf B_{N+2} = \mathsf 0, \\ \mathsf B_k &= \mathsf A \mathsf B_{k+1} - \mathsf B_{k+2} + c_k \mathsf I, \qquad\text{for $N\ge k \ge 1$},\\ \mathsf G(x,y) &= \mathsf B_1 \mathsf F_1(x,y), \end{align} \]
</p>
<p> where \(\mathsf B_k\) are \(2\times2\) matrices. The divided difference \(\Delta[g](x,y)\) is given by the second element of \(\mathsf G(x,y)\).</p>
<p>The same basic recursion applies to the more general case, </p><p class="formulaDsp">
\[ g(x) = \sum_{k=0}^N c_k \sin\bigl( (k+k_0)x + x_0\bigr). \]
</p>
<p> For example, the sum area arising in the computation of geodesic areas, </p><p class="formulaDsp">
\[ \sum_{k=0}^N c_k\cos\bigl((2k+1)\sigma) \]
</p>
<p> is given by \(x=2\sigma\), \(x_0=\frac12\pi\), \(k_0=\frac12\). Again we consider the matrix sum, </p><p class="formulaDsp">
\[ \mathsf G(x,y) = \begin{bmatrix} (g(x) + g(y)) / 2\\ (g(x) - g(y)) / (x - y) \end{bmatrix} = \sum_{k=0}^N c_k \mathsf F_k(x,y), \]
</p>
<p> where </p><p class="formulaDsp">
\[ \mathsf F_k(x,y)= \begin{bmatrix} \cos\bigl( (k+k_0)d \bigr) \sin \bigl( (k+k_0)p + x_0 \bigr)\\ {\displaystyle\frac{\sin\bigl( (k+k_0)d \bigr)}d} \cos \bigl( (k+k_0)p + x_0 \bigr) \end{bmatrix}. \]
</p>
<p> The recursion for \(\mathsf B_k\) is identical to the previous case; in particular, the expression for \(\mathsf A(x,y)\) remains unchanged. The result for the sum is </p><p class="formulaDsp">
\[ \mathsf G(x,y) = (c_0\mathsf I - \mathsf B_2) \mathsf F_0(x,y) + \mathsf B_1 \mathsf F_1(x,y). \]
</p>
<center> Back to <a class="el" href="jacobi.html">Jacobi's conformal projection</a>. Forward to <a class="el" href="greatellipse.html">Great Ellipses</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  <a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/>
</a> 1.8.13
</small></address>
</body>
</html>
|