GeographicLib  1.38
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Rhumb.hpp
Go to the documentation of this file.
1 /**
2  * \file Rhumb.hpp
3  * \brief Header for GeographicLib::Rhumb and GeographicLib::RhumbLine classes
4  *
5  * Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
6  * the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9 
10 #if !defined(GEOGRAPHICLIB_RHUMB_HPP)
11 #define GEOGRAPHICLIB_RHUMB_HPP 1
12 
15 
16 namespace GeographicLib {
17 
18  class RhumbLine;
19 
20  /**
21  * \brief Solve of the direct and inverse rhumb problems.
22  *
23  * The path of constant azimuth between two points on a ellipsoid at (\e
24  * lat1, \e lon1) and (\e lat2, \e lon2) is called the rhumb line (also
25  * called the loxodrome). Its length is \e s12 and its azimuth is \e azi12
26  * and \e azi2. (The azimuth is the heading measured clockwise from north.)
27  *
28  * Given \e lat1, \e lon1, \e azi12, and \e s12, we can determine \e lat2,
29  * and \e lon2. This is the \e direct rhumb problem and its solution is
30  * given by the function Rhumb::Direct.
31  *
32  * Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi12
33  * and \e s12. This is the \e inverse rhumb problem, whose solution is
34  * given by Rhumb::Inverse. This finds the shortest such rhumb line, i.e.,
35  * the one that wraps no more than half way around the earth .
36  *
37  * Note that rhumb lines may be appreciably longer (up to 50%) than the
38  * corresponding Geodesic. For example the distance between London Heathrow
39  * and Tokyo Narita via the rhumb line is 11400 km which is 18% longer than
40  * the geodesic distance 9600 km.
41  *
42  * For more information on rhumb lines see \ref rhumb.
43  *
44  * Example of use:
45  * \include example-Rhumb.cpp
46  **********************************************************************/
47 
49  private:
50  typedef Math::real real;
51  friend class RhumbLine;
52  Ellipsoid _ell;
53  bool _exact;
54  static const int tm_maxord = GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER;
55  static inline real overflow() {
56  // Overflow value s.t. atan(overflow_) = pi/2
57  static const real
58  overflow = 1 / Math::sq(std::numeric_limits<real>::epsilon());
59  return overflow;
60  }
61  static inline real tano(real x) {
62  using std::abs; using std::tan;
63  return
64  2 * abs(x) == Math::pi() ? (x < 0 ? - overflow() : overflow()) :
65  tan(x);
66  }
67  static inline real gd(real x)
68  { using std::atan; using std::sinh; return atan(sinh(x)); }
69 
70  // Use divided differences to determine (mu2 - mu1) / (psi2 - psi1)
71  // accurately
72  //
73  // Definition: Df(x,y,d) = (f(x) - f(y)) / (x - y)
74  // See:
75  // W. M. Kahan and R. J. Fateman,
76  // Symbolic computation of divided differences,
77  // SIGSAM Bull. 33(3), 7-28 (1999)
78  // http://dx.doi.org/10.1145/334714.334716
79  // http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf
80 
81  static inline real Dtan(real x, real y) {
82  real d = x - y, tx = tano(x), ty = tano(y), txy = tx * ty;
83  return d ? (2 * txy > -1 ? (1 + txy) * tano(d) : tx - ty) / d :
84  1 + txy;
85  }
86  static inline real Datan(real x, real y) {
87  using std::atan;
88  real d = x - y, xy = x * y;
89  return d ? (2 * xy > -1 ? atan( d / (1 + xy) ) : atan(x) - atan(y)) / d :
90  1 / (1 + xy);
91  }
92  static inline real Dsin(real x, real y) {
93  using std::sin; using std::cos;
94  real d = (x - y)/2;
95  return cos((x + y)/2) * (d ? sin(d) / d : 1);
96  }
97  static inline real Dsinh(real x, real y) {
98  using std::sinh; using std::cosh;
99  real d = (x - y)/2;
100  return cosh((x + y)/2) * (d ? sinh(d) / d : 1);
101  }
102  static inline real Dasinh(real x, real y) {
103  real d = x - y,
104  hx = Math::hypot(real(1), x), hy = Math::hypot(real(1), y);
105  return d ? Math::asinh(x*y > 0 ? d * (x + y) / (x*hy + y*hx) :
106  x*hy - y*hx) / d :
107  1 / hx;
108  }
109  static inline real Dgd(real x, real y) {
110  using std::sinh;
111  return Datan(sinh(x), sinh(y)) * Dsinh(x, y);
112  }
113  static inline real Dgdinv(real x, real y) {
114  return Dasinh(tano(x), tano(y)) * Dtan(x, y);
115  }
116  // Copied from LambertConformalConic...
117  // e * atanh(e * x) = log( ((1 + e*x)/(1 - e*x))^(e/2) ) if f >= 0
118  // - sqrt(-e2) * atan( sqrt(-e2) * x) if f < 0
119  inline real eatanhe(real x) const {
120  using std::atan;
121  return _ell._f >= 0 ? _ell._e * Math::atanh(_ell._e * x) :
122  - _ell._e * atan(_ell._e * x);
123  }
124  // Copied from LambertConformalConic...
125  // Deatanhe(x,y) = eatanhe((x-y)/(1-e^2*x*y))/(x-y)
126  inline real Deatanhe(real x, real y) const {
127  real t = x - y, d = 1 - _ell._e2 * x * y;
128  return t ? eatanhe(t / d) / t : _ell._e2 / d;
129  }
130  // (E(x) - E(y)) / (x - y) -- E = incomplete elliptic integral of 2nd kind
131  real DE(real x, real y) const;
132  // (mux - muy) / (phix - phiy) using elliptic integrals
133  real DRectifying(real latx, real laty) const;
134  // (psix - psiy) / (phix - phiy)
135  real DIsometric(real latx, real laty) const;
136 
137  // (sum(c[j]*sin(2*j*x),j=1..n) - sum(c[j]*sin(2*j*x),j=1..n)) / (x - y)
138  static real SinSeries(real x, real y, const real c[], int n);
139  // (mux - muy) / (chix - chiy) using Krueger's series
140  real DConformalToRectifying(real chix, real chiy) const;
141  // (chix - chiy) / (mux - muy) using Krueger's series
142  real DRectifyingToConformal(real mux, real muy) const;
143 
144  // (mux - muy) / (psix - psiy)
145  real DIsometricToRectifying(real psix, real psiy) const;
146  // (psix - psiy) / (mux - muy)
147  real DRectifyingToIsometric(real mux, real muy) const;
148 
149  public:
150 
151  /**
152  * Constructor for a ellipsoid with
153  *
154  * @param[in] a equatorial radius (meters).
155  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
156  * Negative \e f gives a prolate ellipsoid. If \e f &gt; 1, set
157  * flattening to 1/\e f.
158  * @param[in] exact if true (the default) use an addition theorem for
159  * elliptic integrals to compute divided differences; otherwise use
160  * series expansion (accurate for |<i>f</i>| < 0.01).
161  * @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
162  * positive.
163  *
164  * See \ref rhumb, for a detailed description of the \e exact parameter.
165  **********************************************************************/
166  Rhumb(real a, real f, bool exact = true) : _ell(a, f), _exact(exact) {}
167 
168  /**
169  * Solve the direct rhumb problem.
170  *
171  * @param[in] lat1 latitude of point 1 (degrees).
172  * @param[in] lon1 longitude of point 1 (degrees).
173  * @param[in] azi12 azimuth of the rhumb line (degrees).
174  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
175  * negative.
176  * @param[out] lat2 latitude of point 2 (degrees).
177  * @param[out] lon2 longitude of point 2 (degrees).
178  *
179  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
180  * azi1 should be in the range [&minus;540&deg;, 540&deg;). The values of
181  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
182  * 180&deg;).
183  *
184  * If point 1 is a pole, the cosine of its latitude is taken to be
185  * 1/&epsilon;<sup>2</sup> (where &epsilon; is 2<sup>-52</sup>). This
186  * position, which is extremely close to the actual pole, allows the
187  * calculation to be carried out in finite terms. If \e s12 is large
188  * enough that the rhumb line crosses a pole, the longitude of point 2
189  * is indeterminate (a NaN is returned for \e lon2).
190  **********************************************************************/
191  void Direct(real lat1, real lon1, real azi12, real s12,
192  real& lat2, real& lon2) const;
193 
194  /**
195  * Solve the inverse rhumb problem.
196  *
197  * @param[in] lat1 latitude of point 1 (degrees).
198  * @param[in] lon1 longitude of point 1 (degrees).
199  * @param[in] lat2 latitude of point 2 (degrees).
200  * @param[in] lon2 longitude of point 2 (degrees).
201  * @param[out] s12 rhumb distance between point 1 and point 2 (meters).
202  * @param[out] azi12 azimuth of the rhumb line (degrees).
203  *
204  * The shortest rhumb line is found. \e lat1 and \e lat2 should be in the
205  * range [&minus;90&deg;, 90&deg;]; \e lon1 and \e lon2 should be in the
206  * range [&minus;540&deg;, 540&deg;). The value of \e azi12 returned is in
207  * the range [&minus;180&deg;, 180&deg;).
208  *
209  * If either point is a pole, the cosine of its latitude is taken to be
210  * 1/&epsilon;<sup>2</sup> (where &epsilon; is 2<sup>-52</sup>). This
211  * position, which is extremely close to the actual pole, allows the
212  * calculation to be carried out in finite terms.
213  **********************************************************************/
214  void Inverse(real lat1, real lon1, real lat2, real lon2,
215  real& s12, real& azi12) const;
216 
217  /**
218  * Set up to compute several points on a single rhumb line.
219  *
220  * @param[in] lat1 latitude of point 1 (degrees).
221  * @param[in] lon1 longitude of point 1 (degrees).
222  * @param[in] azi12 azimuth of the rhumb line (degrees).
223  * @return a RhumbLine object.
224  *
225  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
226  * azi12 should be in the range [&minus;540&deg;, 540&deg;).
227  *
228  * If point 1 is a pole, the cosine of its latitude is taken to be
229  * 1/&epsilon;<sup>2</sup> (where &epsilon; is 2<sup>-52</sup>). This
230  * position, which is extremely close to the actual pole, allows the
231  * calculation to be carried out in finite terms.
232  **********************************************************************/
233  RhumbLine Line(real lat1, real lon1, real azi12) const;
234 
235  /** \name Inspector functions.
236  **********************************************************************/
237  ///@{
238 
239  /**
240  * @return \e a the equatorial radius of the ellipsoid (meters). This is
241  * the value used in the constructor.
242  **********************************************************************/
243  Math::real MajorRadius() const { return _ell.MajorRadius(); }
244 
245  /**
246  * @return \e f the flattening of the ellipsoid. This is the
247  * value used in the constructor.
248  **********************************************************************/
249  Math::real Flattening() const { return _ell.Flattening(); }
250 
251  /**
252  * A global instantiation of Rhumb with the parameters for the WGS84
253  * ellipsoid.
254  **********************************************************************/
255  static const Rhumb& WGS84();
256  };
257 
258  /**
259  * \brief Find a sequence of points on a single rhumb line.
260  *
261  * RhumbLine facilitates the determination of a series of points on a single
262  * rhumb line. The starting point (\e lat1, \e lon1) and the azimuth \e
263  * azi12 are specified in the call to Rhumb::Line which returns a RhumbLine
264  * object. RhumbLine.Position returns the location of point 2 a distance \e
265  * s12 along the rhumb line.
266 
267  * There is no public constructor for this class. (Use Rhumb::Line to create
268  * an instance.) The Rhumb object used to create a RhumbLine must stay in
269  * scope as long as the RhumbLine.
270  *
271  * Example of use:
272  * \include example-RhumbLine.cpp
273  **********************************************************************/
274 
276  private:
277  typedef Math::real real;
278  friend class Rhumb;
279  const Rhumb& _rh;
280  bool _exact;
281  real _lat1, _lon1, _azi12, _salp, _calp, _mu1, _psi1, _r1;
282  RhumbLine& operator=(const RhumbLine&); // copy assignment not allowed
283  RhumbLine(const Rhumb& rh, real lat1, real lon1, real azi12,
284  bool exact);
285  public:
286  /**
287  * Compute the position of point 2 which is a distance \e s12 (meters) from
288  * point 1.
289  *
290  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
291  * negative.
292  * @param[out] lat2 latitude of point 2 (degrees).
293  * @param[out] lon2 longitude of point 2 (degrees).
294  *
295  * The values of \e lon2 and \e azi2 returned are in the range
296  * [&minus;180&deg;, 180&deg;).
297  *
298  * If \e s12 is large enough that the rhumb line crosses a pole, the
299  * longitude of point 2 is indeterminate (a NaN is returned for \e lon2).
300  **********************************************************************/
301  void Position(real s12, real& lat2, real& lon2) const;
302 
303  /** \name Inspector functions
304  **********************************************************************/
305  ///@{
306 
307  /**
308  * @return \e lat1 the latitude of point 1 (degrees).
309  **********************************************************************/
310  Math::real Latitude() const { return _lat1; }
311 
312  /**
313  * @return \e lon1 the longitude of point 1 (degrees).
314  **********************************************************************/
315  Math::real Longitude() const { return _lon1; }
316 
317  /**
318  * @return \e azi12 the azimuth of the rhumb line (degrees).
319  **********************************************************************/
320  Math::real Azimuth() const { return _azi12; }
321 
322  /**
323  * @return \e a the equatorial radius of the ellipsoid (meters). This is
324  * the value inherited from the Rhumb object used in the constructor.
325  **********************************************************************/
326  Math::real MajorRadius() const { return _rh.MajorRadius(); }
327 
328  /**
329  * @return \e f the flattening of the ellipsoid. This is the value
330  * inherited from the Rhumb object used in the constructor.
331  **********************************************************************/
332  Math::real Flattening() const { return _rh.Flattening(); }
333  };
334 
335 } // namespace GeographicLib
336 
337 #endif // GEOGRAPHICLIB_RHUMB_HPP