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