GeographicLib 2.2
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 (2014-2023) <charles@karney.com> and licensed
6 * under the MIT/X11 License. For more information, see
7 * https://geographiclib.sourceforge.io/
8 **********************************************************************/
9
10#if !defined(GEOGRAPHICLIB_RHUMB_HPP)
11#define GEOGRAPHICLIB_RHUMB_HPP 1
12
15#include <vector>
16
17#if !defined(GEOGRAPHICLIB_RHUMBAREA_ORDER)
18/**
19 * The order of the series approximation used in rhumb area calculations.
20 * GEOGRAPHICLIB_RHUMBAREA_ORDER can be set to one of [4, 5, 6, 7, 8].
21 **********************************************************************/
22# define GEOGRAPHICLIB_RHUMBAREA_ORDER \
23 (GEOGRAPHICLIB_PRECISION == 2 ? 6 : \
24 (GEOGRAPHICLIB_PRECISION == 1 ? 4 : 8))
25#endif
26
27#if defined(_MSC_VER)
28// Squelch warnings about dll vs vector
29# pragma warning (push)
30# pragma warning (disable: 4251)
31#endif
32
33namespace GeographicLib {
34
35 class RhumbLine;
36
37 /**
38 * \brief Solve of the direct and inverse rhumb problems.
39 *
40 * The path of constant azimuth between two points on an ellipsoid at (\e
41 * lat1, \e lon1) and (\e lat2, \e lon2) is called the rhumb line (also
42 * called the loxodrome). Its length is \e s12 and its azimuth is \e azi12.
43 * (The azimuth is the heading measured clockwise from north.)
44 *
45 * Given \e lat1, \e lon1, \e azi12, and \e s12, we can determine \e lat2,
46 * and \e lon2. This is the \e direct rhumb problem and its solution is
47 * given by the function Rhumb::Direct.
48 *
49 * Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi12
50 * and \e s12. This is the \e inverse rhumb problem, whose solution is given
51 * by Rhumb::Inverse. This finds the shortest such rhumb line, i.e., the one
52 * that wraps no more than half way around the earth. If the end points are
53 * on opposite meridians, there are two shortest rhumb lines and the
54 * east-going one is chosen.
55 *
56 * These routines also optionally calculate the area under the rhumb line, \e
57 * S12. This is the area, measured counter-clockwise, of the rhumb line
58 * quadrilateral with corners (<i>lat1</i>,<i>lon1</i>), (0,<i>lon1</i>),
59 * (0,<i>lon2</i>), and (<i>lat2</i>,<i>lon2</i>).
60 *
61 * Note that rhumb lines may be appreciably longer (up to 50%) than the
62 * corresponding Geodesic. For example the distance between London Heathrow
63 * and Tokyo Narita via the rhumb line is 11400 km which is 18% longer than
64 * the geodesic distance 9600 km.
65 *
66 * This implementation is described in
67 * - C. F. F. Karney,<br>
68 * <a href="https://arxiv.org/abs/2303.03219">The area of rhumb
69 * polygons</a>,<br>
70 * Technical Report, SRI International, March 2023.<br>
71 * <a href="https://arxiv.org/abs/2303.03219">arxiv:2303.03219</a>
72 * .
73 * For more information on rhumb lines see \ref rhumb.
74 *
75 * Example of use:
76 * \include example-Rhumb.cpp
77 **********************************************************************/
78
80 private:
81 typedef Math::real real;
82 friend class RhumbLine;
83 template<class T> friend class PolygonAreaT;
84 DAuxLatitude _aux;
85 bool _exact;
86 real _a, _f, _n, _rm, _c2;
87 int _lL; // N.B. names of the form _[A-Z].* are reserved in C++
88 std::vector<real> _pP; // The Fourier coefficients P_l
89 static const int Lmax_ = GEOGRAPHICLIB_RHUMBAREA_ORDER;
90 void AreaCoeffs();
91 class qIntegrand {
92 const AuxLatitude& _aux;
93 public:
94 qIntegrand(const AuxLatitude& aux);
95 real operator()(real chi) const;
96 };
97
98 real MeanSinXi(const AuxAngle& chix, const AuxAngle& chiy) const;
99
100 // The following two functions (with lots of ignored arguments) mimic the
101 // interface to the corresponding Geodesic function. These are needed by
102 // PolygonAreaT.
103 void GenDirect(real lat1, real lon1, real azi12,
104 bool, real s12, unsigned outmask,
105 real& lat2, real& lon2, real&, real&, real&, real&, real&,
106 real& S12) const {
107 GenDirect(lat1, lon1, azi12, s12, outmask, lat2, lon2, S12);
108 }
109 void GenInverse(real lat1, real lon1, real lat2, real lon2,
110 unsigned outmask, real& s12, real& azi12,
111 real&, real& , real& , real& , real& S12) const {
112 GenInverse(lat1, lon1, lat2, lon2, outmask, s12, azi12, S12);
113 }
114
115 public:
116 /**
117 * Bit masks for what calculations to do. They specify which results to
118 * return in the general routines Rhumb::GenDirect and Rhumb::GenInverse
119 * routines. RhumbLine::mask is a duplication of this enum.
120 **********************************************************************/
121 enum mask {
122 /**
123 * No output.
124 * @hideinitializer
125 **********************************************************************/
126 NONE = 0U,
127 /**
128 * Calculate latitude \e lat2.
129 * @hideinitializer
130 **********************************************************************/
131 LATITUDE = 1U<<7,
132 /**
133 * Calculate longitude \e lon2.
134 * @hideinitializer
135 **********************************************************************/
136 LONGITUDE = 1U<<8,
137 /**
138 * Calculate azimuth \e azi12.
139 * @hideinitializer
140 **********************************************************************/
141 AZIMUTH = 1U<<9,
142 /**
143 * Calculate distance \e s12.
144 * @hideinitializer
145 **********************************************************************/
146 DISTANCE = 1U<<10,
147 /**
148 * Calculate area \e S12.
149 * @hideinitializer
150 **********************************************************************/
151 AREA = 1U<<14,
152 /**
153 * Unroll \e lon2 in the direct calculation.
154 * @hideinitializer
155 **********************************************************************/
156 LONG_UNROLL = 1U<<15,
157 /**
158 * Calculate everything. (LONG_UNROLL is not included in this mask.)
159 * @hideinitializer
160 **********************************************************************/
161 ALL = 0x7F80U,
162 };
163
164 /**
165 * Constructor for an ellipsoid with
166 *
167 * @param[in] a equatorial radius (meters).
168 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
169 * Negative \e f gives a prolate ellipsoid.
170 * @param[in] exact if true use the exact expressions for the auxiliary
171 * latitudes; otherwise use series expansion (accurate for |<i>f</i>| <
172 * 0.01) [default false].
173 * @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
174 * positive.
175 **********************************************************************/
176 Rhumb(real a, real f, bool exact = false);
177
178 /**
179 * Solve the direct rhumb problem returning also the area.
180 *
181 * @param[in] lat1 latitude of point 1 (degrees).
182 * @param[in] lon1 longitude of point 1 (degrees).
183 * @param[in] azi12 azimuth of the rhumb line (degrees).
184 * @param[in] s12 distance between point 1 and point 2 (meters); it can be
185 * negative.
186 * @param[out] lat2 latitude of point 2 (degrees).
187 * @param[out] lon2 longitude of point 2 (degrees).
188 * @param[out] S12 area under the rhumb line (meters<sup>2</sup>).
189 *
190 * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The value of
191 * \e lon2 returned is in the range [&minus;180&deg;, 180&deg;].
192 *
193 * If point 1 is a pole, the cosine of its latitude is taken to be
194 * 1/&epsilon;<sup>2</sup> (where &epsilon; is 2<sup>-52</sup>). This
195 * position, which is extremely close to the actual pole, allows the
196 * calculation to be carried out in finite terms. If \e s12 is large
197 * enough that the rhumb line crosses a pole, the longitude of point 2
198 * is indeterminate (a NaN is returned for \e lon2 and \e S12).
199 **********************************************************************/
200 void Direct(real lat1, real lon1, real azi12, real s12,
201 real& lat2, real& lon2, real& S12) const {
202 GenDirect(lat1, lon1, azi12, s12,
203 LATITUDE | LONGITUDE | AREA, lat2, lon2, S12);
204 }
205
206 /**
207 * Solve the direct rhumb problem without the area.
208 **********************************************************************/
209 void Direct(real lat1, real lon1, real azi12, real s12,
210 real& lat2, real& lon2) const {
211 real t;
212 GenDirect(lat1, lon1, azi12, s12, LATITUDE | LONGITUDE, lat2, lon2, t);
213 }
214
215 /**
216 * The general direct rhumb problem. Rhumb::Direct is defined in terms
217 * of this function.
218 *
219 * @param[in] lat1 latitude of point 1 (degrees).
220 * @param[in] lon1 longitude of point 1 (degrees).
221 * @param[in] azi12 azimuth of the rhumb line (degrees).
222 * @param[in] s12 distance between point 1 and point 2 (meters); it can be
223 * negative.
224 * @param[in] outmask a bitor'ed combination of Rhumb::mask values
225 * specifying which of the following parameters should be set.
226 * @param[out] lat2 latitude of point 2 (degrees).
227 * @param[out] lon2 longitude of point 2 (degrees).
228 * @param[out] S12 area under the rhumb line (meters<sup>2</sup>).
229 *
230 * The Rhumb::mask values possible for \e outmask are
231 * - \e outmask |= Rhumb::LATITUDE for the latitude \e lat2;
232 * - \e outmask |= Rhumb::LONGITUDE for the latitude \e lon2;
233 * - \e outmask |= Rhumb::AREA for the area \e S12;
234 * - \e outmask |= Rhumb::ALL for all of the above;
235 * - \e outmask |= Rhumb::LONG_UNROLL to unroll \e lon2 instead of wrapping
236 * it into the range [&minus;180&deg;, 180&deg;].
237 * .
238 * With the Rhumb::LONG_UNROLL bit set, the quantity \e lon2 &minus;
239 * \e lon1 indicates how many times and in what sense the rhumb line
240 * encircles the ellipsoid.
241 **********************************************************************/
242 void GenDirect(real lat1, real lon1, real azi12, real s12,
243 unsigned outmask, real& lat2, real& lon2, real& S12) const;
244
245 /**
246 * Solve the inverse rhumb problem returning also the area.
247 *
248 * @param[in] lat1 latitude of point 1 (degrees).
249 * @param[in] lon1 longitude of point 1 (degrees).
250 * @param[in] lat2 latitude of point 2 (degrees).
251 * @param[in] lon2 longitude of point 2 (degrees).
252 * @param[out] s12 rhumb distance between point 1 and point 2 (meters).
253 * @param[out] azi12 azimuth of the rhumb line (degrees).
254 * @param[out] S12 area under the rhumb line (meters<sup>2</sup>).
255 *
256 * The shortest rhumb line is found. If the end points are on opposite
257 * meridians, there are two shortest rhumb lines and the east-going one is
258 * chosen. \e lat1 and \e lat2 should be in the range [&minus;90&deg;,
259 * 90&deg;]. The value of \e azi12 returned is in the range
260 * [&minus;180&deg;, 180&deg;].
261 *
262 * If either point is a pole, the cosine of its latitude is taken to be
263 * 1/&epsilon;<sup>2</sup> (where &epsilon; is 2<sup>-52</sup>). This
264 * position, which is extremely close to the actual pole, allows the
265 * calculation to be carried out in finite terms.
266 **********************************************************************/
267 void Inverse(real lat1, real lon1, real lat2, real lon2,
268 real& s12, real& azi12, real& S12) const {
269 GenInverse(lat1, lon1, lat2, lon2,
270 DISTANCE | AZIMUTH | AREA, s12, azi12, S12);
271 }
272
273 /**
274 * Solve the inverse rhumb problem without the area.
275 **********************************************************************/
276 void Inverse(real lat1, real lon1, real lat2, real lon2,
277 real& s12, real& azi12) const {
278 real t;
279 GenInverse(lat1, lon1, lat2, lon2, DISTANCE | AZIMUTH, s12, azi12, t);
280 }
281
282 /**
283 * The general inverse rhumb problem. Rhumb::Inverse is defined in terms
284 * of this function.
285 *
286 * @param[in] lat1 latitude of point 1 (degrees).
287 * @param[in] lon1 longitude of point 1 (degrees).
288 * @param[in] lat2 latitude of point 2 (degrees).
289 * @param[in] lon2 longitude of point 2 (degrees).
290 * @param[in] outmask a bitor'ed combination of Rhumb::mask values
291 * specifying which of the following parameters should be set.
292 * @param[out] s12 rhumb distance between point 1 and point 2 (meters).
293 * @param[out] azi12 azimuth of the rhumb line (degrees).
294 * @param[out] S12 area under the rhumb line (meters<sup>2</sup>).
295 *
296 * The Rhumb::mask values possible for \e outmask are
297 * - \e outmask |= Rhumb::DISTANCE for the latitude \e s12;
298 * - \e outmask |= Rhumb::AZIMUTH for the latitude \e azi12;
299 * - \e outmask |= Rhumb::AREA for the area \e S12;
300 * - \e outmask |= Rhumb::ALL for all of the above;
301 **********************************************************************/
302 void GenInverse(real lat1, real lon1, real lat2, real lon2,
303 unsigned outmask,
304 real& s12, real& azi12, real& S12) const;
305
306 /**
307 * Set up to compute several points on a single rhumb line.
308 *
309 * @param[in] lat1 latitude of point 1 (degrees).
310 * @param[in] lon1 longitude of point 1 (degrees).
311 * @param[in] azi12 azimuth of the rhumb line (degrees).
312 * @return a RhumbLine object.
313 *
314 * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
315 *
316 * If point 1 is a pole, the cosine of its latitude is taken to be
317 * 1/&epsilon;<sup>2</sup> (where &epsilon; is 2<sup>-52</sup>). This
318 * position, which is extremely close to the actual pole, allows the
319 * calculation to be carried out in finite terms.
320 **********************************************************************/
321 RhumbLine Line(real lat1, real lon1, real azi12) const;
322
323 /** \name Inspector functions.
324 **********************************************************************/
325 ///@{
326
327 /**
328 * @return \e a the equatorial radius of the ellipsoid (meters). This is
329 * the value used in the constructor.
330 **********************************************************************/
331 Math::real EquatorialRadius() const { return _a; }
332
333 /**
334 * @return \e f the flattening of the ellipsoid. This is the
335 * value used in the constructor.
336 **********************************************************************/
337 Math::real Flattening() const { return _f; }
338
339 /**
340 * @return total area of ellipsoid in meters<sup>2</sup>. The area of a
341 * polygon encircling a pole can be found by adding
342 * Geodesic::EllipsoidArea()/2 to the sum of \e S12 for each side of the
343 * polygon.
344 **********************************************************************/
346 // _c2 contains a Math::degrees() factor, so 4*pi -> 2*Math::td.
347 return 2 * real(Math::td) * _c2;
348 }
349 ///@}
350
351 /**
352 * A global instantiation of Rhumb with the parameters for the WGS84
353 * ellipsoid.
354 **********************************************************************/
355 static const Rhumb& WGS84();
356 };
357
358 /**
359 * \brief Find a sequence of points on a single rhumb line.
360 *
361 * RhumbLine facilitates the determination of a series of points on a single
362 * rhumb line. The starting point (\e lat1, \e lon1) and the azimuth \e
363 * azi12 are specified in the call to Rhumb::Line which returns a RhumbLine
364 * object. RhumbLine.Position returns the location of point 2 (and,
365 * optionally, the corresponding area, \e S12) a distance \e s12 along the
366 * rhumb line.
367 *
368 * There is no public constructor for this class. (Use Rhumb::Line to create
369 * an instance.) The Rhumb object used to create a RhumbLine must stay in
370 * scope as long as the RhumbLine.
371 *
372 * Example of use:
373 * \include example-RhumbLine.cpp
374 **********************************************************************/
375
377 private:
378 typedef Math::real real;
379 friend class Rhumb;
380 const Rhumb& _rh;
381 real _lat1, _lon1, _azi12, _salp, _calp, _mu1, _psi1;
382 AuxAngle _phi1, _chi1;
383 // copy assignment not allowed
384 RhumbLine& operator=(const RhumbLine&) = delete;
385 RhumbLine(const Rhumb& rh, real lat1, real lon1, real azi12);
386
387 public:
388 /**
389 * Construction is via default copy constructor.
390 **********************************************************************/
391 RhumbLine(const RhumbLine&) = default;
392 /**
393 * This is a duplication of Rhumb::mask.
394 **********************************************************************/
395 enum mask {
396 /**
397 * No output.
398 * @hideinitializer
399 **********************************************************************/
400 NONE = Rhumb::NONE,
401 /**
402 * Calculate latitude \e lat2.
403 * @hideinitializer
404 **********************************************************************/
405 LATITUDE = Rhumb::LATITUDE,
406 /**
407 * Calculate longitude \e lon2.
408 * @hideinitializer
409 **********************************************************************/
410 LONGITUDE = Rhumb::LONGITUDE,
411 /**
412 * Calculate azimuth \e azi12.
413 * @hideinitializer
414 **********************************************************************/
415 AZIMUTH = Rhumb::AZIMUTH,
416 /**
417 * Calculate distance \e s12.
418 * @hideinitializer
419 **********************************************************************/
420 DISTANCE = Rhumb::DISTANCE,
421 /**
422 * Calculate area \e S12.
423 * @hideinitializer
424 **********************************************************************/
425 AREA = Rhumb::AREA,
426 /**
427 * Unroll \e lon2 in the direct calculation.
428 * @hideinitializer
429 **********************************************************************/
430 LONG_UNROLL = Rhumb::LONG_UNROLL,
431 /**
432 * Calculate everything. (LONG_UNROLL is not included in this mask.)
433 * @hideinitializer
434 **********************************************************************/
435 ALL = Rhumb::ALL,
436 };
437
438 /**
439 * Compute the position of point 2 which is a distance \e s12 (meters) from
440 * point 1. The area is also computed.
441 *
442 * @param[in] s12 distance between point 1 and point 2 (meters); it can be
443 * negative.
444 * @param[out] lat2 latitude of point 2 (degrees).
445 * @param[out] lon2 longitude of point 2 (degrees).
446 * @param[out] S12 area under the rhumb line (meters<sup>2</sup>).
447 *
448 * The value of \e lon2 returned is in the range [&minus;180&deg;,
449 * 180&deg;].
450 *
451 * If \e s12 is large enough that the rhumb line crosses a pole, the
452 * longitude of point 2 is indeterminate (a NaN is returned for \e lon2 and
453 * \e S12).
454 **********************************************************************/
455 void Position(real s12, real& lat2, real& lon2, real& S12) const {
456 GenPosition(s12, LATITUDE | LONGITUDE | AREA, lat2, lon2, S12);
457 }
458
459 /**
460 * Compute the position of point 2 which is a distance \e s12 (meters) from
461 * point 1. The area is not computed.
462 **********************************************************************/
463 void Position(real s12, real& lat2, real& lon2) const {
464 real t;
465 GenPosition(s12, LATITUDE | LONGITUDE, lat2, lon2, t);
466 }
467
468 /**
469 * The general position routine. RhumbLine::Position is defined in term so
470 * this function.
471 *
472 * @param[in] s12 distance between point 1 and point 2 (meters); it can be
473 * negative.
474 * @param[in] outmask a bitor'ed combination of RhumbLine::mask values
475 * specifying which of the following parameters should be set.
476 * @param[out] lat2 latitude of point 2 (degrees).
477 * @param[out] lon2 longitude of point 2 (degrees).
478 * @param[out] S12 area under the rhumb line (meters<sup>2</sup>).
479 *
480 * The RhumbLine::mask values possible for \e outmask are
481 * - \e outmask |= RhumbLine::LATITUDE for the latitude \e lat2;
482 * - \e outmask |= RhumbLine::LONGITUDE for the latitude \e lon2;
483 * - \e outmask |= RhumbLine::AREA for the area \e S12;
484 * - \e outmask |= RhumbLine::ALL for all of the above;
485 * - \e outmask |= RhumbLine::LONG_UNROLL to unroll \e lon2 instead of
486 * wrapping it into the range [&minus;180&deg;, 180&deg;].
487 * .
488 * With the RhumbLine::LONG_UNROLL bit set, the quantity \e lon2 &minus; \e
489 * lon1 indicates how many times and in what sense the rhumb line encircles
490 * the ellipsoid.
491 *
492 * If \e s12 is large enough that the rhumb line crosses a pole, the
493 * longitude of point 2 is indeterminate (a NaN is returned for \e lon2 and
494 * \e S12).
495 **********************************************************************/
496 void GenPosition(real s12, unsigned outmask,
497 real& lat2, real& lon2, real& S12) const;
498
499 /** \name Inspector functions
500 **********************************************************************/
501 ///@{
502
503 /**
504 * @return \e lat1 the latitude of point 1 (degrees).
505 **********************************************************************/
506 Math::real Latitude() const { return _lat1; }
507
508 /**
509 * @return \e lon1 the longitude of point 1 (degrees).
510 **********************************************************************/
511 Math::real Longitude() const { return _lon1; }
512
513 /**
514 * @return \e azi12 the azimuth of the rhumb line (degrees).
515 **********************************************************************/
516 Math::real Azimuth() const { return _azi12; }
517
518 /**
519 * @return \e a the equatorial radius of the ellipsoid (meters). This is
520 * the value inherited from the Rhumb object used in the constructor.
521 **********************************************************************/
523
524 /**
525 * @return \e f the flattening of the ellipsoid. This is the value
526 * inherited from the Rhumb object used in the constructor.
527 **********************************************************************/
528 Math::real Flattening() const { return _rh.Flattening(); }
529 };
530
531} // namespace GeographicLib
532
533#endif // GEOGRAPHICLIB_RHUMB_HPP
Header for GeographicLib::Constants class.
#define GEOGRAPHICLIB_EXPORT
Definition Constants.hpp:67
Header for the GeographicLib::DAuxLatitude class.
GeographicLib::Math::real real
Definition GeodSolve.cpp:31
#define GEOGRAPHICLIB_RHUMBAREA_ORDER
Definition Rhumb.hpp:22
An accurate representation of angles.
Definition AuxAngle.hpp:43
Conversions between auxiliary latitudes.
Divided differences of auxiliary latitudes.
Find a sequence of points on a single rhumb line.
Definition Rhumb.hpp:376
void Position(real s12, real &lat2, real &lon2) const
Definition Rhumb.hpp:463
void Position(real s12, real &lat2, real &lon2, real &S12) const
Definition Rhumb.hpp:455
RhumbLine(const RhumbLine &)=default
Math::real EquatorialRadius() const
Definition Rhumb.hpp:522
Math::real Latitude() const
Definition Rhumb.hpp:506
Math::real Azimuth() const
Definition Rhumb.hpp:516
Math::real Flattening() const
Definition Rhumb.hpp:528
Math::real Longitude() const
Definition Rhumb.hpp:511
Solve of the direct and inverse rhumb problems.
Definition Rhumb.hpp:79
Math::real Flattening() const
Definition Rhumb.hpp:337
void Direct(real lat1, real lon1, real azi12, real s12, real &lat2, real &lon2, real &S12) const
Definition Rhumb.hpp:200
void Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi12, real &S12) const
Definition Rhumb.hpp:267
void Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi12) const
Definition Rhumb.hpp:276
void Direct(real lat1, real lon1, real azi12, real s12, real &lat2, real &lon2) const
Definition Rhumb.hpp:209
Math::real EquatorialRadius() const
Definition Rhumb.hpp:331
Math::real EllipsoidArea() const
Definition Rhumb.hpp:345
Namespace for GeographicLib.