GeographicLib  1.41
Geodesic.hpp
Go to the documentation of this file.
1 /**
2  * \file Geodesic.hpp
3  * \brief Header for GeographicLib::Geodesic class
4  *
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9
10 #if !defined(GEOGRAPHICLIB_GEODESIC_HPP)
11 #define GEOGRAPHICLIB_GEODESIC_HPP 1
12
14
15 #if !defined(GEOGRAPHICLIB_GEODESIC_ORDER)
16 /**
17  * The order of the expansions used by Geodesic.
18  * GEOGRAPHICLIB_GEODESIC_ORDER can be set to any integer in [0, 8].
19  **********************************************************************/
20 # define GEOGRAPHICLIB_GEODESIC_ORDER \
21  (GEOGRAPHICLIB_PRECISION == 2 ? 6 : \
22  (GEOGRAPHICLIB_PRECISION == 1 ? 3 : \
23  (GEOGRAPHICLIB_PRECISION == 3 ? 7 : 8)))
24 #endif
25
26 namespace GeographicLib {
27
28  class GeodesicLine;
29
30  /**
31  * \brief %Geodesic calculations
32  *
33  * The shortest path between two points on a ellipsoid at (\e lat1, \e lon1)
34  * and (\e lat2, \e lon2) is called the geodesic. Its length is \e s12 and
35  * the geodesic from point 1 to point 2 has azimuths \e azi1 and \e azi2 at
36  * the two end points. (The azimuth is the heading measured clockwise from
37  * north. \e azi2 is the "forward" azimuth, i.e., the heading that takes you
38  * beyond point 2 not back to point 1.) In the figure below, latitude if
39  * labeled &phi;, longitude &lambda; (with &lambda;<sub>12</sub> =
40  * &lambda;<sub>2</sub> &minus; &lambda;<sub>1</sub>), and azimuth &alpha;.
41  *
42  * <img src="http://upload.wikimedia.org/wikipedia/commons/c/cb/Geodesic_problem_on_an_ellipsoid.svg" width=250 alt="spheroidal triangle">
43  *
44  * Given \e lat1, \e lon1, \e azi1, and \e s12, we can determine \e lat2, \e
45  * lon2, and \e azi2. This is the \e direct geodesic problem and its
46  * solution is given by the function Geodesic::Direct. (If \e s12 is
47  * sufficiently large that the geodesic wraps more than halfway around the
48  * earth, there will be another geodesic between the points with a smaller \e
49  * s12.)
50  *
51  * Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi1, \e
52  * azi2, and \e s12. This is the \e inverse geodesic problem, whose solution
53  * is given by Geodesic::Inverse. Usually, the solution to the inverse
54  * problem is unique. In cases where there are multiple solutions (all with
55  * the same \e s12, of course), all the solutions can be easily generated
56  * once a particular solution is provided.
57  *
58  * The standard way of specifying the direct problem is the specify the
59  * distance \e s12 to the second point. However it is sometimes useful
60  * instead to specify the arc length \e a12 (in degrees) on the auxiliary
61  * sphere. This is a mathematical construct used in solving the geodesic
62  * problems. The solution of the direct problem in this form is provided by
63  * Geodesic::ArcDirect. An arc length in excess of 180&deg; indicates that
64  * the geodesic is not a shortest path. In addition, the arc length between
65  * an equatorial crossing and the next extremum of latitude for a geodesic is
66  * 90&deg;.
67  *
68  * This class can also calculate several other quantities related to
69  * geodesics. These are:
70  * - <i>reduced length</i>. If we fix the first point and increase \e azi1
71  * by \e dazi1 (radians), the second point is displaced \e m12 \e dazi1 in
72  * the direction \e azi2 + 90&deg;. The quantity \e m12 is called
73  * the "reduced length" and is symmetric under interchange of the two
74  * points. On a curved surface the reduced length obeys a symmetry
75  * relation, \e m12 + \e m21 = 0. On a flat surface, we have \e m12 = \e
76  * s12. The ratio <i>s12</i>/\e m12 gives the azimuthal scale for an
77  * azimuthal equidistant projection.
78  * - <i>geodesic scale</i>. Consider a reference geodesic and a second
79  * geodesic parallel to this one at point 1 and separated by a small
80  * distance \e dt. The separation of the two geodesics at point 2 is \e
81  * M12 \e dt where \e M12 is called the "geodesic scale". \e M21 is
82  * defined similarly (with the geodesics being parallel at point 2). On a
83  * flat surface, we have \e M12 = \e M21 = 1. The quantity 1/\e M12 gives
84  * the scale of the Cassini-Soldner projection.
85  * - <i>area</i>. The area between the geodesic from point 1 to point 2 and
86  * the equation is represented by \e S12; it is the area, measured
87  * counter-clockwise, of the geodesic quadrilateral with corners
88  * (<i>lat1</i>,<i>lon1</i>), (0,<i>lon1</i>), (0,<i>lon2</i>), and
89  * (<i>lat2</i>,<i>lon2</i>). It can be used to compute the area of any
90  * simple geodesic polygon.
91  *
92  * Overloaded versions of Geodesic::Direct, Geodesic::ArcDirect, and
93  * Geodesic::Inverse allow these quantities to be returned. In addition
94  * there are general functions Geodesic::GenDirect, and Geodesic::GenInverse
95  * which allow an arbitrary set of results to be computed. The quantities \e
96  * m12, \e M12, \e M21 which all specify the behavior of nearby geodesics
97  * obey addition rules. If points 1, 2, and 3 all lie on a single geodesic,
98  * then the following rules hold:
99  * - \e s13 = \e s12 + \e s23
100  * - \e a13 = \e a12 + \e a23
101  * - \e S13 = \e S12 + \e S23
102  * - \e m13 = \e m12 \e M23 + \e m23 \e M21
103  * - \e M13 = \e M12 \e M23 &minus; (1 &minus; \e M12 \e M21) \e m23 / \e m12
104  * - \e M31 = \e M32 \e M21 &minus; (1 &minus; \e M23 \e M32) \e m12 / \e m23
105  *
106  * Additional functionality is provided by the GeodesicLine class, which
107  * allows a sequence of points along a geodesic to be computed.
108  *
109  * The shortest distance returned by the solution of the inverse problem is
110  * (obviously) uniquely defined. However, in a few special cases there are
111  * multiple azimuths which yield the same shortest distance. Here is a
112  * catalog of those cases:
113  * - \e lat1 = &minus;\e lat2 (with neither point at a pole). If \e azi1 =
114  * \e azi2, the geodesic is unique. Otherwise there are two geodesics and
115  * the second one is obtained by setting [\e azi1, \e azi2] = [\e azi2, \e
116  * azi1], [\e M12, \e M21] = [\e M21, \e M12], \e S12 = &minus;\e S12.
117  * (This occurs when the longitude difference is near &plusmn;180&deg; for
118  * oblate ellipsoids.)
119  * - \e lon2 = \e lon1 &plusmn; 180&deg; (with neither point at a pole). If
120  * \e azi1 = 0&deg; or &plusmn;180&deg;, the geodesic is unique. Otherwise
121  * there are two geodesics and the second one is obtained by setting [\e
122  * azi1, \e azi2] = [&minus;\e azi1, &minus;\e azi2], \e S12 = &minus;\e
123  * S12. (This occurs when \e lat2 is near &minus;\e lat1 for prolate
124  * ellipsoids.)
125  * - Points 1 and 2 at opposite poles. There are infinitely many geodesics
126  * which can be generated by setting [\e azi1, \e azi2] = [\e azi1, \e
127  * azi2] + [\e d, &minus;\e d], for arbitrary \e d. (For spheres, this
128  * prescription applies when points 1 and 2 are antipodal.)
129  * - s12 = 0 (coincident points). There are infinitely many geodesics which
130  * can be generated by setting [\e azi1, \e azi2] = [\e azi1, \e azi2] +
131  * [\e d, \e d], for arbitrary \e d.
132  *
133  * The calculations are accurate to better than 15 nm (15 nanometers) for the
134  * WGS84 ellipsoid. See Sec. 9 of
135  * <a href="http://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a> for
136  * details. The algorithms used by this class are based on series expansions
137  * using the flattening \e f as a small parameter. These are only accurate
138  * for |<i>f</i>| &lt; 0.02; however reasonably accurate results will be
139  * obtained for |<i>f</i>| &lt; 0.2. Here is a table of the approximate
140  * maximum error (expressed as a distance) for an ellipsoid with the same
141  * major radius as the WGS84 ellipsoid and different values of the
142  * flattening.<pre>
143  * |f| error
144  * 0.01 25 nm
145  * 0.02 30 nm
146  * 0.05 10 um
147  * 0.1 1.5 mm
148  * 0.2 300 mm
149  * </pre>
150  * For very eccentric ellipsoids, use GeodesicExact instead.
151  *
152  * The algorithms are described in
153  * - C. F. F. Karney,
154  * <a href="https://dx.doi.org/10.1007/s00190-012-0578-z">
155  * Algorithms for geodesics</a>,
156  * J. Geodesy <b>87</b>, 43--55 (2013);
157  * DOI: <a href="https://dx.doi.org/10.1007/s00190-012-0578-z">
158  * 10.1007/s00190-012-0578-z</a>;
161  * .
163  *
164  * Example of use:
165  * \include example-Geodesic.cpp
166  *
167  * <a href="GeodSolve.1.html">GeodSolve</a> is a command-line utility
169  **********************************************************************/
170
172  private:
173  typedef Math::real real;
174  friend class GeodesicLine;
175  static const int nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
176  static const int nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
177  static const int nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER;
178  static const int nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
179  static const int nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
180  static const int nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
181  static const int nA3x_ = nA3_;
182  static const int nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
183  static const int nC3x_ = (nC3_ * (nC3_ - 1)) / 2;
184  static const int nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER;
185  static const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
186  static const unsigned maxit1_ = 20;
187  unsigned maxit2_;
188  real tiny_, tol0_, tol1_, tol2_, tolb_, xthresh_;
189
190  enum captype {
191  CAP_NONE = 0U,
192  CAP_C1 = 1U<<0,
193  CAP_C1p = 1U<<1,
194  CAP_C2 = 1U<<2,
195  CAP_C3 = 1U<<3,
196  CAP_C4 = 1U<<4,
197  CAP_ALL = 0x1FU,
199  OUT_ALL = 0x7F80U,
200  OUT_MASK = 0xFF80U, // Includes LONG_NOWRAP
201  };
202
203  static real SinCosSeries(bool sinp,
204  real sinx, real cosx, const real c[], int n);
205  static inline real AngRound(real x) {
206  // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
207  // for reals = 0.7 pm on the earth if x is an angle in degrees. (This
208  // is about 1000 times more resolution than we get with angles around 90
209  // degrees.) We use this to avoid having to deal with near singular
210  // cases when x is non-zero but tiny (e.g., 1.0e-200).
211  using std::abs;
212  const real z = 1/real(16);
213  GEOGRAPHICLIB_VOLATILE real y = abs(x);
214  // The compiler mustn't "simplify" z - (z - y) to y
215  y = y < z ? z - (z - y) : y;
216  return x < 0 ? -y : y;
217  }
218  static inline void SinCosNorm(real& sinx, real& cosx) {
219  real r = Math::hypot(sinx, cosx);
220  sinx /= r;
221  cosx /= r;
222  }
223  static real Astroid(real x, real y);
224
225  real _a, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2;
226  real _A3x[nA3x_], _C3x[nC3x_], _C4x[nC4x_];
227
228  void Lengths(real eps, real sig12,
229  real ssig1, real csig1, real dn1,
230  real ssig2, real csig2, real dn2,
231  real cbet1, real cbet2,
232  real& s12s, real& m12a, real& m0,
233  bool scalep, real& M12, real& M21,
234  real C1a[], real C2a[]) const;
235  real InverseStart(real sbet1, real cbet1, real dn1,
236  real sbet2, real cbet2, real dn2,
237  real lam12,
238  real& salp1, real& calp1,
239  real& salp2, real& calp2, real& dnm,
240  real C1a[], real C2a[]) const;
241  real Lambda12(real sbet1, real cbet1, real dn1,
242  real sbet2, real cbet2, real dn2,
243  real salp1, real calp1,
244  real& salp2, real& calp2, real& sig12,
245  real& ssig1, real& csig1, real& ssig2, real& csig2,
246  real& eps, real& domg12, bool diffp, real& dlam12,
247  real C1a[], real C2a[], real C3a[])
248  const;
249
250  // These are Maxima generated functions to provide series approximations to
251  // the integrals for the ellipsoidal geodesic.
252  static real A1m1f(real eps);
253  static void C1f(real eps, real c[]);
254  static void C1pf(real eps, real c[]);
255  static real A2m1f(real eps);
256  static void C2f(real eps, real c[]);
257
258  void A3coeff();
259  real A3f(real eps) const;
260  void C3coeff();
261  void C3f(real eps, real c[]) const;
262  void C4coeff();
263  void C4f(real k2, real c[]) const;
264
265  public:
266
267  /**
268  * Bit masks for what calculations to do. These masks do double duty.
269  * They signify to the GeodesicLine::GeodesicLine constructor and to
270  * Geodesic::Line what capabilities should be included in the GeodesicLine
271  * object. They also specify which results to return in the general
272  * routines Geodesic::GenDirect and Geodesic::GenInverse routines.
273  * GeodesicLine::mask is a duplication of this enum.
274  **********************************************************************/
276  /**
277  * No capabilities, no output.
278  * @hideinitializer
279  **********************************************************************/
280  NONE = 0U,
281  /**
282  * Calculate latitude \e lat2. (It's not necessary to include this as a
283  * capability to GeodesicLine because this is included by default.)
284  * @hideinitializer
285  **********************************************************************/
286  LATITUDE = 1U<<7 | CAP_NONE,
287  /**
288  * Calculate longitude \e lon2.
289  * @hideinitializer
290  **********************************************************************/
291  LONGITUDE = 1U<<8 | CAP_C3,
292  /**
293  * Calculate azimuths \e azi1 and \e azi2. (It's not necessary to
294  * include this as a capability to GeodesicLine because this is included
295  * by default.)
296  * @hideinitializer
297  **********************************************************************/
298  AZIMUTH = 1U<<9 | CAP_NONE,
299  /**
300  * Calculate distance \e s12.
301  * @hideinitializer
302  **********************************************************************/
303  DISTANCE = 1U<<10 | CAP_C1,
304  /**
305  * Allow distance \e s12 to be used as input in the direct geodesic
306  * problem.
307  * @hideinitializer
308  **********************************************************************/
309  DISTANCE_IN = 1U<<11 | CAP_C1 | CAP_C1p,
310  /**
311  * Calculate reduced length \e m12.
312  * @hideinitializer
313  **********************************************************************/
314  REDUCEDLENGTH = 1U<<12 | CAP_C1 | CAP_C2,
315  /**
316  * Calculate geodesic scales \e M12 and \e M21.
317  * @hideinitializer
318  **********************************************************************/
319  GEODESICSCALE = 1U<<13 | CAP_C1 | CAP_C2,
320  /**
321  * Calculate area \e S12.
322  * @hideinitializer
323  **********************************************************************/
324  AREA = 1U<<14 | CAP_C4,
325  /**
326  * Do not wrap the \e lon2 in the direct calculation.
327  * @hideinitializer
328  **********************************************************************/
329  LONG_NOWRAP = 1U<<15,
330  /**
331  * All capabilities, calculate everything. (LONG_NOWRAP is not
332  * included in this mask.)
333  * @hideinitializer
334  **********************************************************************/
335  ALL = OUT_ALL| CAP_ALL,
336  };
337
338  /** \name Constructor
339  **********************************************************************/
340  ///@{
341  /**
342  * Constructor for a ellipsoid with
343  *
344  * @param[in] a equatorial radius (meters).
345  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
346  * Negative \e f gives a prolate ellipsoid. If \e f &gt; 1, set
347  * flattening to 1/\e f.
348  * @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
349  * positive.
350  **********************************************************************/
351  Geodesic(real a, real f);
352  ///@}
353
354  /** \name Direct geodesic problem specified in terms of distance.
355  **********************************************************************/
356  ///@{
357  /**
358  * Solve the direct geodesic problem where the length of the geodesic
359  * is specified in terms of distance.
360  *
361  * @param[in] lat1 latitude of point 1 (degrees).
362  * @param[in] lon1 longitude of point 1 (degrees).
363  * @param[in] azi1 azimuth at point 1 (degrees).
364  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
365  * negative.
366  * @param[out] lat2 latitude of point 2 (degrees).
367  * @param[out] lon2 longitude of point 2 (degrees).
368  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
369  * @param[out] m12 reduced length of geodesic (meters).
370  * @param[out] M12 geodesic scale of point 2 relative to point 1
371  * (dimensionless).
372  * @param[out] M21 geodesic scale of point 1 relative to point 2
373  * (dimensionless).
374  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
375  * @return \e a12 arc length of between point 1 and point 2 (degrees).
376  *
377  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
378  * azi1 should be in the range [&minus;540&deg;, 540&deg;). The values of
379  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
380  * 180&deg;).
381  *
382  * If either point is at a pole, the azimuth is defined by keeping the
383  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
384  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
385  * 180&deg; signifies a geodesic which is not a shortest path. (For a
386  * prolate ellipsoid, an additional condition is necessary for a shortest
387  * path: the longitudinal extent must not exceed of 180&deg;.)
388  *
389  * The following functions are overloaded versions of Geodesic::Direct
390  * which omit some of the output parameters. Note, however, that the arc
391  * length is always computed and returned as the function value.
392  **********************************************************************/
393  Math::real Direct(real lat1, real lon1, real azi1, real s12,
394  real& lat2, real& lon2, real& azi2,
395  real& m12, real& M12, real& M21, real& S12)
396  const {
397  real t;
398  return GenDirect(lat1, lon1, azi1, false, s12,
399  LATITUDE | LONGITUDE | AZIMUTH |
400  REDUCEDLENGTH | GEODESICSCALE | AREA,
401  lat2, lon2, azi2, t, m12, M12, M21, S12);
402  }
403
404  /**
405  * See the documentation for Geodesic::Direct.
406  **********************************************************************/
407  Math::real Direct(real lat1, real lon1, real azi1, real s12,
408  real& lat2, real& lon2)
409  const {
410  real t;
411  return GenDirect(lat1, lon1, azi1, false, s12,
412  LATITUDE | LONGITUDE,
413  lat2, lon2, t, t, t, t, t, t);
414  }
415
416  /**
417  * See the documentation for Geodesic::Direct.
418  **********************************************************************/
419  Math::real Direct(real lat1, real lon1, real azi1, real s12,
420  real& lat2, real& lon2, real& azi2)
421  const {
422  real t;
423  return GenDirect(lat1, lon1, azi1, false, s12,
424  LATITUDE | LONGITUDE | AZIMUTH,
425  lat2, lon2, azi2, t, t, t, t, t);
426  }
427
428  /**
429  * See the documentation for Geodesic::Direct.
430  **********************************************************************/
431  Math::real Direct(real lat1, real lon1, real azi1, real s12,
432  real& lat2, real& lon2, real& azi2, real& m12)
433  const {
434  real t;
435  return GenDirect(lat1, lon1, azi1, false, s12,
436  LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH,
437  lat2, lon2, azi2, t, m12, t, t, t);
438  }
439
440  /**
441  * See the documentation for Geodesic::Direct.
442  **********************************************************************/
443  Math::real Direct(real lat1, real lon1, real azi1, real s12,
444  real& lat2, real& lon2, real& azi2,
445  real& M12, real& M21)
446  const {
447  real t;
448  return GenDirect(lat1, lon1, azi1, false, s12,
449  LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE,
450  lat2, lon2, azi2, t, t, M12, M21, t);
451  }
452
453  /**
454  * See the documentation for Geodesic::Direct.
455  **********************************************************************/
456  Math::real Direct(real lat1, real lon1, real azi1, real s12,
457  real& lat2, real& lon2, real& azi2,
458  real& m12, real& M12, real& M21)
459  const {
460  real t;
461  return GenDirect(lat1, lon1, azi1, false, s12,
462  LATITUDE | LONGITUDE | AZIMUTH |
463  REDUCEDLENGTH | GEODESICSCALE,
464  lat2, lon2, azi2, t, m12, M12, M21, t);
465  }
466  ///@}
467
468  /** \name Direct geodesic problem specified in terms of arc length.
469  **********************************************************************/
470  ///@{
471  /**
472  * Solve the direct geodesic problem where the length of the geodesic
473  * is specified in terms of arc length.
474  *
475  * @param[in] lat1 latitude of point 1 (degrees).
476  * @param[in] lon1 longitude of point 1 (degrees).
477  * @param[in] azi1 azimuth at point 1 (degrees).
478  * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
479  * be negative.
480  * @param[out] lat2 latitude of point 2 (degrees).
481  * @param[out] lon2 longitude of point 2 (degrees).
482  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
483  * @param[out] s12 distance between point 1 and point 2 (meters).
484  * @param[out] m12 reduced length of geodesic (meters).
485  * @param[out] M12 geodesic scale of point 2 relative to point 1
486  * (dimensionless).
487  * @param[out] M21 geodesic scale of point 1 relative to point 2
488  * (dimensionless).
489  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
490  *
491  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
492  * azi1 should be in the range [&minus;540&deg;, 540&deg;). The values of
493  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
494  * 180&deg;).
495  *
496  * If either point is at a pole, the azimuth is defined by keeping the
497  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
498  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
499  * 180&deg; signifies a geodesic which is not a shortest path. (For a
500  * prolate ellipsoid, an additional condition is necessary for a shortest
501  * path: the longitudinal extent must not exceed of 180&deg;.)
502  *
503  * The following functions are overloaded versions of Geodesic::Direct
504  * which omit some of the output parameters.
505  **********************************************************************/
506  void ArcDirect(real lat1, real lon1, real azi1, real a12,
507  real& lat2, real& lon2, real& azi2, real& s12,
508  real& m12, real& M12, real& M21, real& S12)
509  const {
510  GenDirect(lat1, lon1, azi1, true, a12,
511  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
512  REDUCEDLENGTH | GEODESICSCALE | AREA,
513  lat2, lon2, azi2, s12, m12, M12, M21, S12);
514  }
515
516  /**
517  * See the documentation for Geodesic::ArcDirect.
518  **********************************************************************/
519  void ArcDirect(real lat1, real lon1, real azi1, real a12,
520  real& lat2, real& lon2) const {
521  real t;
522  GenDirect(lat1, lon1, azi1, true, a12,
523  LATITUDE | LONGITUDE,
524  lat2, lon2, t, t, t, t, t, t);
525  }
526
527  /**
528  * See the documentation for Geodesic::ArcDirect.
529  **********************************************************************/
530  void ArcDirect(real lat1, real lon1, real azi1, real a12,
531  real& lat2, real& lon2, real& azi2) const {
532  real t;
533  GenDirect(lat1, lon1, azi1, true, a12,
534  LATITUDE | LONGITUDE | AZIMUTH,
535  lat2, lon2, azi2, t, t, t, t, t);
536  }
537
538  /**
539  * See the documentation for Geodesic::ArcDirect.
540  **********************************************************************/
541  void ArcDirect(real lat1, real lon1, real azi1, real a12,
542  real& lat2, real& lon2, real& azi2, real& s12)
543  const {
544  real t;
545  GenDirect(lat1, lon1, azi1, true, a12,
546  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
547  lat2, lon2, azi2, s12, t, t, t, t);
548  }
549
550  /**
551  * See the documentation for Geodesic::ArcDirect.
552  **********************************************************************/
553  void ArcDirect(real lat1, real lon1, real azi1, real a12,
554  real& lat2, real& lon2, real& azi2,
555  real& s12, real& m12) const {
556  real t;
557  GenDirect(lat1, lon1, azi1, true, a12,
558  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
559  REDUCEDLENGTH,
560  lat2, lon2, azi2, s12, m12, t, t, t);
561  }
562
563  /**
564  * See the documentation for Geodesic::ArcDirect.
565  **********************************************************************/
566  void ArcDirect(real lat1, real lon1, real azi1, real a12,
567  real& lat2, real& lon2, real& azi2, real& s12,
568  real& M12, real& M21) const {
569  real t;
570  GenDirect(lat1, lon1, azi1, true, a12,
571  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
572  GEODESICSCALE,
573  lat2, lon2, azi2, s12, t, M12, M21, t);
574  }
575
576  /**
577  * See the documentation for Geodesic::ArcDirect.
578  **********************************************************************/
579  void ArcDirect(real lat1, real lon1, real azi1, real a12,
580  real& lat2, real& lon2, real& azi2, real& s12,
581  real& m12, real& M12, real& M21) const {
582  real t;
583  GenDirect(lat1, lon1, azi1, true, a12,
584  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
585  REDUCEDLENGTH | GEODESICSCALE,
586  lat2, lon2, azi2, s12, m12, M12, M21, t);
587  }
588  ///@}
589
590  /** \name General version of the direct geodesic solution.
591  **********************************************************************/
592  ///@{
593
594  /**
595  * The general direct geodesic problem. Geodesic::Direct and
596  * Geodesic::ArcDirect are defined in terms of this function.
597  *
598  * @param[in] lat1 latitude of point 1 (degrees).
599  * @param[in] lon1 longitude of point 1 (degrees).
600  * @param[in] azi1 azimuth at point 1 (degrees).
601  * @param[in] arcmode boolean flag determining the meaning of the \e
602  * s12_a12.
603  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
604  * point 1 and point 2 (meters); otherwise it is the arc length between
605  * point 1 and point 2 (degrees); it can be negative.
607  * specifying which of the following parameters should be set.
608  * @param[out] lat2 latitude of point 2 (degrees).
609  * @param[out] lon2 longitude of point 2 (degrees).
610  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
611  * @param[out] s12 distance between point 1 and point 2 (meters).
612  * @param[out] m12 reduced length of geodesic (meters).
613  * @param[out] M12 geodesic scale of point 2 relative to point 1
614  * (dimensionless).
615  * @param[out] M21 geodesic scale of point 1 relative to point 2
616  * (dimensionless).
617  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
618  * @return \e a12 arc length of between point 1 and point 2 (degrees).
619  *
621  * - \e outmask |= Geodesic::LATITUDE for the latitude \e lat2;
622  * - \e outmask |= Geodesic::LONGITUDE for the latitude \e lon2;
623  * - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2;
624  * - \e outmask |= Geodesic::DISTANCE for the distance \e s12;
625  * - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e
626  * m12;
627  * - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e
628  * M12 and \e M21;
629  * - \e outmask |= Geodesic::AREA for the area \e S12;
630  * - \e outmask |= Geodesic::ALL for all of the above;
631  * - \e outmask |= Geodesic::LONG_NOWRAP stops the returned value of \e
632  * lon2 being wrapped into the range [&minus;180&deg;, 180&deg;).
633  * .
634  * The function value \e a12 is always computed and returned and this
635  * equals \e s12_a12 is \e arcmode is true. If \e outmask includes
636  * Geodesic::DISTANCE and \e arcmode is false, then \e s12 = \e s12_a12.
637  * It is not necessary to include Geodesic::DISTANCE_IN in \e outmask; this
638  * is automatically included is \e arcmode is false.
639  *
640  * With the LONG_NOWRAP bit set, the quantity \e lon2 &minus; \e lon1
641  * indicates how many times the geodesic wrapped around the ellipsoid.
642  * Because \e lon2 might be outside the normal allowed range for
643  * longitudes, [&minus;540&deg;, 540&deg;), be sure to normalize it with
644  * Math::AngNormalize2 before using it in other GeographicLib calls.
645  **********************************************************************/
646  Math::real GenDirect(real lat1, real lon1, real azi1,
647  bool arcmode, real s12_a12, unsigned outmask,
648  real& lat2, real& lon2, real& azi2,
649  real& s12, real& m12, real& M12, real& M21,
650  real& S12) const;
651  ///@}
652
653  /** \name Inverse geodesic problem.
654  **********************************************************************/
655  ///@{
656  /**
657  * Solve the inverse geodesic problem.
658  *
659  * @param[in] lat1 latitude of point 1 (degrees).
660  * @param[in] lon1 longitude of point 1 (degrees).
661  * @param[in] lat2 latitude of point 2 (degrees).
662  * @param[in] lon2 longitude of point 2 (degrees).
663  * @param[out] s12 distance between point 1 and point 2 (meters).
664  * @param[out] azi1 azimuth at point 1 (degrees).
665  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
666  * @param[out] m12 reduced length of geodesic (meters).
667  * @param[out] M12 geodesic scale of point 2 relative to point 1
668  * (dimensionless).
669  * @param[out] M21 geodesic scale of point 1 relative to point 2
670  * (dimensionless).
671  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
672  * @return \e a12 arc length of between point 1 and point 2 (degrees).
673  *
674  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e
675  * lon1 and \e lon2 should be in the range [&minus;540&deg;, 540&deg;).
676  * The values of \e azi1 and \e azi2 returned are in the range
677  * [&minus;180&deg;, 180&deg;).
678  *
679  * If either point is at a pole, the azimuth is defined by keeping the
680  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
681  * and taking the limit &epsilon; &rarr; 0+.
682  *
683  * The solution to the inverse problem is found using Newton's method. If
684  * this fails to converge (this is very unlikely in geodetic applications
685  * but does occur for very eccentric ellipsoids), then the bisection method
686  * is used to refine the solution.
687  *
688  * The following functions are overloaded versions of Geodesic::Inverse
689  * which omit some of the output parameters. Note, however, that the arc
690  * length is always computed and returned as the function value.
691  **********************************************************************/
692  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
693  real& s12, real& azi1, real& azi2, real& m12,
694  real& M12, real& M21, real& S12) const {
695  return GenInverse(lat1, lon1, lat2, lon2,
696  DISTANCE | AZIMUTH |
697  REDUCEDLENGTH | GEODESICSCALE | AREA,
698  s12, azi1, azi2, m12, M12, M21, S12);
699  }
700
701  /**
702  * See the documentation for Geodesic::Inverse.
703  **********************************************************************/
704  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
705  real& s12) const {
706  real t;
707  return GenInverse(lat1, lon1, lat2, lon2,
708  DISTANCE,
709  s12, t, t, t, t, t, t);
710  }
711
712  /**
713  * See the documentation for Geodesic::Inverse.
714  **********************************************************************/
715  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
716  real& azi1, real& azi2) const {
717  real t;
718  return GenInverse(lat1, lon1, lat2, lon2,
719  AZIMUTH,
720  t, azi1, azi2, t, t, t, t);
721  }
722
723  /**
724  * See the documentation for Geodesic::Inverse.
725  **********************************************************************/
726  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
727  real& s12, real& azi1, real& azi2)
728  const {
729  real t;
730  return GenInverse(lat1, lon1, lat2, lon2,
731  DISTANCE | AZIMUTH,
732  s12, azi1, azi2, t, t, t, t);
733  }
734
735  /**
736  * See the documentation for Geodesic::Inverse.
737  **********************************************************************/
738  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
739  real& s12, real& azi1, real& azi2, real& m12)
740  const {
741  real t;
742  return GenInverse(lat1, lon1, lat2, lon2,
743  DISTANCE | AZIMUTH | REDUCEDLENGTH,
744  s12, azi1, azi2, m12, t, t, t);
745  }
746
747  /**
748  * See the documentation for Geodesic::Inverse.
749  **********************************************************************/
750  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
751  real& s12, real& azi1, real& azi2,
752  real& M12, real& M21) const {
753  real t;
754  return GenInverse(lat1, lon1, lat2, lon2,
755  DISTANCE | AZIMUTH | GEODESICSCALE,
756  s12, azi1, azi2, t, M12, M21, t);
757  }
758
759  /**
760  * See the documentation for Geodesic::Inverse.
761  **********************************************************************/
762  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
763  real& s12, real& azi1, real& azi2, real& m12,
764  real& M12, real& M21) const {
765  real t;
766  return GenInverse(lat1, lon1, lat2, lon2,
767  DISTANCE | AZIMUTH |
768  REDUCEDLENGTH | GEODESICSCALE,
769  s12, azi1, azi2, m12, M12, M21, t);
770  }
771  ///@}
772
773  /** \name General version of inverse geodesic solution.
774  **********************************************************************/
775  ///@{
776  /**
777  * The general inverse geodesic calculation. Geodesic::Inverse is defined
778  * in terms of this function.
779  *
780  * @param[in] lat1 latitude of point 1 (degrees).
781  * @param[in] lon1 longitude of point 1 (degrees).
782  * @param[in] lat2 latitude of point 2 (degrees).
783  * @param[in] lon2 longitude of point 2 (degrees).
785  * specifying which of the following parameters should be set.
786  * @param[out] s12 distance between point 1 and point 2 (meters).
787  * @param[out] azi1 azimuth at point 1 (degrees).
788  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
789  * @param[out] m12 reduced length of geodesic (meters).
790  * @param[out] M12 geodesic scale of point 2 relative to point 1
791  * (dimensionless).
792  * @param[out] M21 geodesic scale of point 1 relative to point 2
793  * (dimensionless).
794  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
795  * @return \e a12 arc length of between point 1 and point 2 (degrees).
796  *
798  * - \e outmask |= Geodesic::DISTANCE for the distance \e s12;
799  * - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2;
800  * - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e
801  * m12;
802  * - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e
803  * M12 and \e M21;
804  * - \e outmask |= Geodesic::AREA for the area \e S12;
805  * - \e outmask |= Geodesic::ALL for all of the above.
806  * .
807  * The arc length is always computed and returned as the function value.
808  **********************************************************************/
809  Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
811  real& s12, real& azi1, real& azi2,
812  real& m12, real& M12, real& M21, real& S12)
813  const;
814  ///@}
815
816  /** \name Interface to GeodesicLine.
817  **********************************************************************/
818  ///@{
819
820  /**
821  * Set up to compute several points on a single geodesic.
822  *
823  * @param[in] lat1 latitude of point 1 (degrees).
824  * @param[in] lon1 longitude of point 1 (degrees).
825  * @param[in] azi1 azimuth at point 1 (degrees).
826  * @param[in] caps bitor'ed combination of Geodesic::mask values
827  * specifying the capabilities the GeodesicLine object should possess,
828  * i.e., which quantities can be returned in calls to
829  * GeodesicLine::Position.
830  * @return a GeodesicLine object.
831  *
832  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
833  * azi1 should be in the range [&minus;540&deg;, 540&deg;).
834  *
835  * The Geodesic::mask values are
836  * - \e caps |= Geodesic::LATITUDE for the latitude \e lat2; this is
838  * - \e caps |= Geodesic::LONGITUDE for the latitude \e lon2;
839  * - \e caps |= Geodesic::AZIMUTH for the azimuth \e azi2; this is
841  * - \e caps |= Geodesic::DISTANCE for the distance \e s12;
842  * - \e caps |= Geodesic::REDUCEDLENGTH for the reduced length \e m12;
843  * - \e caps |= Geodesic::GEODESICSCALE for the geodesic scales \e M12
844  * and \e M21;
845  * - \e caps |= Geodesic::AREA for the area \e S12;
846  * - \e caps |= Geodesic::DISTANCE_IN permits the length of the
847  * geodesic to be given in terms of \e s12; without this capability the
848  * length can only be specified in terms of arc length;
849  * - \e caps |= Geodesic::ALL for all of the above.
850  * .
851  * The default value of \e caps is Geodesic::ALL.
852  *
853  * If the point is at a pole, the azimuth is defined by keeping \e lon1
854  * fixed, writing \e lat1 = &plusmn;(90 &minus; &epsilon;), and taking the
855  * limit &epsilon; &rarr; 0+.
856  **********************************************************************/
857  GeodesicLine Line(real lat1, real lon1, real azi1, unsigned caps = ALL)
858  const;
859
860  ///@}
861
862  /** \name Inspector functions.
863  **********************************************************************/
864  ///@{
865
866  /**
867  * @return \e a the equatorial radius of the ellipsoid (meters). This is
868  * the value used in the constructor.
869  **********************************************************************/
870  Math::real MajorRadius() const { return _a; }
871
872  /**
873  * @return \e f the flattening of the ellipsoid. This is the
874  * value used in the constructor.
875  **********************************************************************/
876  Math::real Flattening() const { return _f; }
877
878  /// \cond SKIP
879  /**
880  * <b>DEPRECATED</b>
881  * @return \e r the inverse flattening of the ellipsoid.
882  **********************************************************************/
883  Math::real InverseFlattening() const { return 1/_f; }
884  /// \endcond
885
886  /**
887  * @return total area of ellipsoid in meters<sup>2</sup>. The area of a
888  * polygon encircling a pole can be found by adding
889  * Geodesic::EllipsoidArea()/2 to the sum of \e S12 for each side of the
890  * polygon.
891  **********************************************************************/
893  { return 4 * Math::pi() * _c2; }
894  ///@}
895
896  /**
897  * A global instantiation of Geodesic with the parameters for the WGS84
898  * ellipsoid.
899  **********************************************************************/
900  static const Geodesic& WGS84();
901
902  };
903
904 } // namespace GeographicLib
905
906 #endif // GEOGRAPHICLIB_GEODESIC_HPP
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &M12, real &M21) const
Definition: Geodesic.hpp:566
static T pi()
Definition: Math.hpp:214
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:71
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &m12) const
Definition: Geodesic.hpp:553
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12, real &M12, real &M21, real &S12) const
Definition: Geodesic.hpp:393
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12) const
Definition: Geodesic.hpp:431
GeographicLib::Math::real real
Definition: GeodSolve.cpp:32
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12, real &M12, real &M21, real &S12) const
Definition: Geodesic.hpp:692
Definition: Geodesic.hpp:870
#define GEOGRAPHICLIB_VOLATILE
Definition: Math.hpp:84
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2) const
Definition: Geodesic.hpp:519
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &m12, real &M12, real &M21) const
Definition: Geodesic.hpp:579
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2) const
Definition: Geodesic.hpp:407
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &M12, real &M21) const
Definition: Geodesic.hpp:750
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12) const
Definition: Geodesic.hpp:738
static T hypot(T x, T y)
Definition: Math.hpp:255
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12, real &M12, real &M21) const
Definition: Geodesic.hpp:456
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2) const
Definition: Geodesic.hpp:726
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12) const
Definition: Geodesic.hpp:704
Math::real EllipsoidArea() const
Definition: Geodesic.hpp:892
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2) const
Definition: Geodesic.hpp:530
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &m12, real &M12, real &M21, real &S12) const
Definition: Geodesic.hpp:506
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12, real &M12, real &M21) const
Definition: Geodesic.hpp:762
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &azi1, real &azi2) const
Definition: Geodesic.hpp:715
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12) const
Definition: Geodesic.hpp:541
#define GEOGRAPHICLIB_GEODESIC_ORDER
Definition: Geodesic.hpp:20
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2) const
Definition: Geodesic.hpp:419
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &M12, real &M21) const
Definition: Geodesic.hpp:443
Math::real Flattening() const
Definition: Geodesic.hpp:876
Geodesic calculations
Definition: Geodesic.hpp:171