- Generated by
1.9.4
|
GeographicLib 2.1.2
|
The Rhumb and RhumbLine classes together with the RhumbSolve 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:
Go to
References:
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
\[ \begin{align} \frac{d\psi}{d\phi} &= \frac{\rho}R, \\ \frac{d\mu}{d\phi} &= \frac{\pi}{2M}\rho, \\ \end{align} \]
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).
Rhumb 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
\[ \tan\alpha_{12} = \frac{\lambda_{12}}{\psi_{12}}, \]
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 shortest rhumb line).
The distance is given by
\[ \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} \]
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
\[ \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} \]
where the last relation is given by the defining equations for \(\psi\) and \(\mu\).
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 Rhumb class, the ellipsoidal generalization is handled by the Ellipsoid class where the auxiliary latitudes are defined.
Here we brief develop the necessary formulas for the Auxiliary latitudes.
Isometric latitude: The equation for \(\psi\) can be integrated to give
\[ \psi = \sinh^{-1}\tan\phi - e\tanh^{-1}(e \sin\phi), \]
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)
\[ \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} \]
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
\[ \begin{align} \phi &= \tan^{-1}\tau,\\ \psi &= \sinh^{-1}\tau'. \end{align} \]
This allows conversions to and from \(\psi\) to be carried out for any value of the flattening \(f\); these conversions are implemented by Ellipsoid::IsometricLatitude and Ellipsoid::InverseIsometricLatitude. (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.)
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.
Parametric latitude: This is given by
\[ \tan\beta = (1-f)\tan\phi, \]
which allows rapid and accurate conversions; these conversions are implemented by Ellipsoid::ParametricLatitude and Ellipsoid::InverseParametricLatitude.
Rectifying latitude: 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\).
In terms of \(\beta\), the rectifying latitude is given by
\[ \mu = \frac{\pi}{2E(ie')} E(\beta,ie'), \]
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 https://dlmf.nist.gov/19.2.ii. These can be evaluated accurately for arbitrary flattening using the method given in https://dlmf.nist.gov/19.36.i. 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 Ellipsoid::RectifyingLatitude and Ellipsoid::InverseRectifyingLatitude.
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 Geodesic class where the coefficients are \(C_{1j}\) and \(C_{1j}'\); see Eqs. (18) and (21) of Karney (2013) and Expansions for geodesics. (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.
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\).
The area between a rhumb line and the equator is given by
\[ S_{12} = \int_{\lambda_1}^{\lambda_2} c^2 \sin\xi \,d\lambda, \]
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
\[ \begin{align} \psi &= m(\lambda-\lambda_0),\\ m &= \frac{\psi_2 - \psi_1}{\lambda_2 - \lambda_1}. \end{align} \]
Performing the integral over \(\lambda\) gives
\[ S_{12} = c^2 \lambda_{12} \left<\sin\xi\right>_{12}, \]
where \(\left<\sin\xi\right>_{12}\) is the mean value of \(\sin\xi\) given by
\[ \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} \]
\(\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.
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
\[ \left<\sin\xi\right>_{12} \rightarrow \frac{dS(\chi_1)/d\chi_1}{\sec\chi_1} =\sin\xi_1, \]
as expected.
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 next section) and use Clenshaw summation to evaluate the sum (see the last section).
Here is the series expansion accurate to 10th order, found by the Maxima script rhumbarea.mac:
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;
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.
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).
Here I give a more systematic treatment using the algebra of divided differences. Many readers will be already using this technique, e.g., when writing, for example,
\[ \frac{\sin(2x) - \sin(2y)}{x-y} = 2\frac{\sin(x-y)}{x-y}\cos(x+y). \]
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 LambertConformalConic and AlbersEqualArea use this technique to improve the accuracy.)
To illustrate the technique, consider the relation between \(\psi\) and \(\chi\),
\[ \psi = \sinh^{-1}\tan\chi. \]
The divided difference operator is defined by
\[ \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} \]
Many of the rules for differentiation apply to divided differences; in particular, we can use the chain rule to write
\[ \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). \]
Kahan and Fateman catalog the divided difference formulas for all the elementary functions. In this case, we need
\[ \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} \]
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.)
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 TransverseMercator 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).
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 https://dlmf.nist.gov/19.11.E2. This can be converted in the following divided difference formula
\[ \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} \]
where the angle \(z\) is given by
\[ \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} \]
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.
The end result is that the Rhumb 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.
The use of Clenshaw summation for summing series of the form,
\[ g(x) = \sum_{k=1}^N c_k \sin kx, \]
is well established. However when computing divided differences, we are interested in evaluating
\[ 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). \]
Clenshaw summation can be used in this case if we simultaneously compute \(g(x)+g(y)\) and perform a matrix summation,
\[ \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), \]
where
\[ \mathsf F_k(x,y)= \begin{bmatrix} \cos kd \sin kp\\ {\displaystyle\frac{\sin kd}d} \cos kp \end{bmatrix}, \]
\(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
\[ \mathsf F_{k+1}(x,y) = \mathsf A(x,y) \mathsf F_k(x,y) - \mathsf F_{k-1}(x,y), \]
where
\[ \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}, \]
and \(\lim_{d\rightarrow0}(\sin d)/d = 1\). The standard Clenshaw algorithm can now be applied to yield
\[ \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} \]
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)\).
The same basic recursion applies to the more general case,
\[ g(x) = \sum_{k=0}^N c_k \sin\bigl( (k+k_0)x + x_0\bigr). \]
For example, the sum area arising in the computation of geodesic areas,
\[ \sum_{k=0}^N c_k\cos\bigl((2k+1)\sigma) \]
is given by \(x=2\sigma\), \(x_0=\frac12\pi\), \(k_0=\frac12\). Again we consider the matrix sum,
\[ \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), \]
where
\[ \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}. \]
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
\[ \mathsf G(x,y) = (c_0\mathsf I - \mathsf B_2) \mathsf F_0(x,y) + \mathsf B_1 \mathsf F_1(x,y). \]
1.9.4