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