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 (this serves to define the course of rhumb lines), the rectifying latitude
(needed for computing distances along rhumb lines), and the parametric latitude
(needed for dealing with rhumb lines that run along a parallel). These are defined in terms of the geographical latitude
by
where is the meridional radius of curvature,
is the radius of a circle of latitude, and
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 on the ellipsoid to a point
. The azimuth
of a rhumb line from
to
is thus given by
where the quadrant of is determined by the signs of the numerator and denominator of the right-hand side,
,
, and, typically,
is reduced to the range
(thus giving the course of the shortest rhumb line).
The distance is given by
where . This relation is indeterminate if
, so we take the limits
and
and apply L'Hôpital's rule to yield
where the last relation is given by the defining equations for and
.
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 can be integrated to give
where is the eccentricity of the ellipsoid. To invert this equation (to give
in terms of
), it is convenient to introduce the variables
and
(
is the conformal latitude) which are related by (Karney, 2011)
The equation for can be inverted to give
in terms of
using Newton's method with
as a starting guess; and, of course,
and
are then given by
This allows conversions to and from to be carried out for any value of the flattening
; these conversions are implemented by Ellipsoid::IsometricLatitude and Ellipsoid::InverseIsometricLatitude. (For prolate ellipsoids,
is negative and
is imaginary, and the equation for
needs to be recast in terms of the tangent function.)
For small values of , Engsager and Poder (2007) express
as a trigonometric series in
. This series can be reverted to give a series for
in terms of
. 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
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 converge faster than the corresponding ones in terms of
.
In terms of , the rectifying latitude is given by
where is the second eccentricity and
and
are the complete and incomplete elliptic integrals of the second kind; see http://dlmf.nist.gov/19.2.ii. These can be evaluated accurately for arbitrary flattening using the method given in http://dlmf.nist.gov/19.36.i. To find
in terms of
, Newton's method can be used (noting that
); for a starting guess use
(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, can be expressed as a series in various ways. The most economical series is in terms of the second flattening
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
given
). These series are used by the Geodesic class where the coefficients are
and
; see Eqs. (18) and (21) of Karney (2013) and Expansions for geodesics. (Make the replacements,
,
,
, 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 ,
,
,
,
to be carried out simply and accurately. The approaches using series apply only if
is small. The others apply for arbitrary values of
.
Despite our ability to compute latitudes accurately, the way that distances enter into the solution involves the ratio ; the numerical calculation of this term is subject to catastrophic round-off errors when
and
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 to be computed accurately (but the formula for
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,
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 and
,
The divided difference operator is defined by
Many of the rules for differentiation apply to divided differences; in particular, we can use the chain rule to write
Kahan and Fateman catalog the divided difference formulas for all the elementary functions. In this case, we need
The crucial point in these formulas is that the right hand sides can be evaluated accurately even if is small. (I've only included the basic results here; Kahan and Fateman supplement these with the derivatives if
vanishes or the direct ratios if
is not small.)
To complete the computation of , we now need
, 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
and
. 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 . If we want to extend the method to arbitrary flattening we need to compute
. The necessary relation is the "addition
theorem" for the incomplete elliptic integral of the second kind given in http://dlmf.nist.gov/19.11.E2. This can be converted in the followinging divided difference formula
where the angle is given by
We also need to apply the divided difference formulas to the conversions from to
and
to
; 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,
is well established. However when computing divided differences, we are interested in evaluating
Clenshaw summation can be used in this case if we simultaneously compute and perform a matrix summation,
where
,
, and, in the limit
,
. The first element of
is the average value of
and the second element, the divided difference, is the average slope.
satisfies the recurrence relation
where
and . The standard Clenshaw algorithm can now be applied to yield
where are
matrices. The divided difference
is given by the second element of
.