GeographicLib  1.46
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-2015) <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 [3, 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] &rarr; [\e
116  * azi2, \e azi1], [\e M12, \e M21] &rarr; [\e M21, \e M12], \e S12 &rarr;
117  * &minus;\e S12. (This occurs when the longitude difference is near
118  * &plusmn;180&deg; for 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] &rarr; [&minus;\e azi1, &minus;\e azi2], \e S12 &rarr;
123  * &minus;\e S12. (This occurs when \e lat2 is near &minus;\e lat1 for
124  * prolate 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] &rarr; [\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] &rarr; [\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  * .
162  * For more information on geodesics see \ref geodesic.
163  *
164  * Example of use:
165  * \include example-Geodesic.cpp
166  *
167  * <a href="GeodSolve.1.html">GeodSolve</a> is a command-line utility
168  * providing access to the functionality of Geodesic and GeodesicLine.
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  // Size for temporary array
187  // nC = max(max(nC1_, nC1p_, nC2_) + 1, max(nC3_, nC4_))
188  static const int nC_ = GEOGRAPHICLIB_GEODESIC_ORDER + 1;
189  static const unsigned maxit1_ = 20;
190  unsigned maxit2_;
191  real tiny_, tol0_, tol1_, tol2_, tolb_, xthresh_;
192
193  enum captype {
194  CAP_NONE = 0U,
195  CAP_C1 = 1U<<0,
196  CAP_C1p = 1U<<1,
197  CAP_C2 = 1U<<2,
198  CAP_C3 = 1U<<3,
199  CAP_C4 = 1U<<4,
200  CAP_ALL = 0x1FU,
201  CAP_MASK = CAP_ALL,
202  OUT_ALL = 0x7F80U,
203  OUT_MASK = 0xFF80U, // Includes LONG_UNROLL
204  };
205
206  static real SinCosSeries(bool sinp,
207  real sinx, real cosx, const real c[], int n);
208  static real Astroid(real x, real y);
209
210  real _a, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2;
211  real _A3x[nA3x_], _C3x[nC3x_], _C4x[nC4x_];
212
213  void Lengths(real eps, real sig12,
214  real ssig1, real csig1, real dn1,
215  real ssig2, real csig2, real dn2,
216  real cbet1, real cbet2, unsigned outmask,
217  real& s12s, real& m12a, real& m0,
218  real& M12, real& M21, real Ca[]) const;
219  real InverseStart(real sbet1, real cbet1, real dn1,
220  real sbet2, real cbet2, real dn2,
221  real lam12, real slam12, real clam12,
222  real& salp1, real& calp1,
223  real& salp2, real& calp2, real& dnm,
224  real Ca[]) const;
225  real Lambda12(real sbet1, real cbet1, real dn1,
226  real sbet2, real cbet2, real dn2,
227  real salp1, real calp1, real slam120, real clam120,
228  real& salp2, real& calp2, real& sig12,
229  real& ssig1, real& csig1, real& ssig2, real& csig2,
230  real& eps, real& somg12, real& comg12,
231  bool diffp, real& dlam12, real Ca[]) const;
232  real GenInverse(real lat1, real lon1, real lat2, real lon2,
233  unsigned outmask, real& s12,
234  real& salp1, real& calp1, real& salp2, real& calp2,
235  real& m12, real& M12, real& M21, real& S12) const;
236
237  // These are Maxima generated functions to provide series approximations to
238  // the integrals for the ellipsoidal geodesic.
239  static real A1m1f(real eps);
240  static void C1f(real eps, real c[]);
241  static void C1pf(real eps, real c[]);
242  static real A2m1f(real eps);
243  static void C2f(real eps, real c[]);
244
245  void A3coeff();
246  real A3f(real eps) const;
247  void C3coeff();
248  void C3f(real eps, real c[]) const;
249  void C4coeff();
250  void C4f(real k2, real c[]) const;
251
252  public:
253
254  /**
255  * Bit masks for what calculations to do. These masks do double duty.
256  * They signify to the GeodesicLine::GeodesicLine constructor and to
257  * Geodesic::Line what capabilities should be included in the GeodesicLine
258  * object. They also specify which results to return in the general
259  * routines Geodesic::GenDirect and Geodesic::GenInverse routines.
260  * GeodesicLine::mask is a duplication of this enum.
261  **********************************************************************/
262  enum mask {
263  /**
264  * No capabilities, no output.
265  * @hideinitializer
266  **********************************************************************/
267  NONE = 0U,
268  /**
269  * Calculate latitude \e lat2. (It's not necessary to include this as a
270  * capability to GeodesicLine because this is included by default.)
271  * @hideinitializer
272  **********************************************************************/
273  LATITUDE = 1U<<7 | CAP_NONE,
274  /**
275  * Calculate longitude \e lon2.
276  * @hideinitializer
277  **********************************************************************/
278  LONGITUDE = 1U<<8 | CAP_C3,
279  /**
280  * Calculate azimuths \e azi1 and \e azi2. (It's not necessary to
281  * include this as a capability to GeodesicLine because this is included
282  * by default.)
283  * @hideinitializer
284  **********************************************************************/
285  AZIMUTH = 1U<<9 | CAP_NONE,
286  /**
287  * Calculate distance \e s12.
288  * @hideinitializer
289  **********************************************************************/
290  DISTANCE = 1U<<10 | CAP_C1,
291  /**
292  * Allow distance \e s12 to be used as input in the direct geodesic
293  * problem.
294  * @hideinitializer
295  **********************************************************************/
296  DISTANCE_IN = 1U<<11 | CAP_C1 | CAP_C1p,
297  /**
298  * Calculate reduced length \e m12.
299  * @hideinitializer
300  **********************************************************************/
301  REDUCEDLENGTH = 1U<<12 | CAP_C1 | CAP_C2,
302  /**
303  * Calculate geodesic scales \e M12 and \e M21.
304  * @hideinitializer
305  **********************************************************************/
306  GEODESICSCALE = 1U<<13 | CAP_C1 | CAP_C2,
307  /**
308  * Calculate area \e S12.
309  * @hideinitializer
310  **********************************************************************/
311  AREA = 1U<<14 | CAP_C4,
312  /**
313  * Unroll \e lon2 in the direct calculation.
314  * @hideinitializer
315  **********************************************************************/
316  LONG_UNROLL = 1U<<15,
317  /**
318  * All capabilities, calculate everything. (LONG_UNROLL is not
319  * included in this mask.)
320  * @hideinitializer
321  **********************************************************************/
322  ALL = OUT_ALL| CAP_ALL,
323  };
324
325  /** \name Constructor
326  **********************************************************************/
327  ///@{
328  /**
329  * Constructor for a ellipsoid with
330  *
331  * @param[in] a equatorial radius (meters).
332  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
333  * Negative \e f gives a prolate ellipsoid.
334  * @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
335  * positive.
336  **********************************************************************/
337  Geodesic(real a, real f);
338  ///@}
339
340  /** \name Direct geodesic problem specified in terms of distance.
341  **********************************************************************/
342  ///@{
343  /**
344  * Solve the direct geodesic problem where the length of the geodesic
345  * is specified in terms of distance.
346  *
347  * @param[in] lat1 latitude of point 1 (degrees).
348  * @param[in] lon1 longitude of point 1 (degrees).
349  * @param[in] azi1 azimuth at point 1 (degrees).
350  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
351  * negative.
352  * @param[out] lat2 latitude of point 2 (degrees).
353  * @param[out] lon2 longitude of point 2 (degrees).
354  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
355  * @param[out] m12 reduced length of geodesic (meters).
356  * @param[out] M12 geodesic scale of point 2 relative to point 1
357  * (dimensionless).
358  * @param[out] M21 geodesic scale of point 1 relative to point 2
359  * (dimensionless).
360  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
361  * @return \e a12 arc length of between point 1 and point 2 (degrees).
362  *
363  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The values of
364  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
365  * 180&deg;).
366  *
367  * If either point is at a pole, the azimuth is defined by keeping the
368  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
369  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
370  * 180&deg; signifies a geodesic which is not a shortest path. (For a
371  * prolate ellipsoid, an additional condition is necessary for a shortest
372  * path: the longitudinal extent must not exceed of 180&deg;.)
373  *
374  * The following functions are overloaded versions of Geodesic::Direct
375  * which omit some of the output parameters. Note, however, that the arc
376  * length is always computed and returned as the function value.
377  **********************************************************************/
378  Math::real Direct(real lat1, real lon1, real azi1, real s12,
379  real& lat2, real& lon2, real& azi2,
380  real& m12, real& M12, real& M21, real& S12)
381  const {
382  real t;
383  return GenDirect(lat1, lon1, azi1, false, s12,
384  LATITUDE | LONGITUDE | AZIMUTH |
385  REDUCEDLENGTH | GEODESICSCALE | AREA,
386  lat2, lon2, azi2, t, m12, M12, M21, S12);
387  }
388
389  /**
390  * See the documentation for Geodesic::Direct.
391  **********************************************************************/
392  Math::real Direct(real lat1, real lon1, real azi1, real s12,
393  real& lat2, real& lon2)
394  const {
395  real t;
396  return GenDirect(lat1, lon1, azi1, false, s12,
397  LATITUDE | LONGITUDE,
398  lat2, lon2, t, t, t, t, t, t);
399  }
400
401  /**
402  * See the documentation for Geodesic::Direct.
403  **********************************************************************/
404  Math::real Direct(real lat1, real lon1, real azi1, real s12,
405  real& lat2, real& lon2, real& azi2)
406  const {
407  real t;
408  return GenDirect(lat1, lon1, azi1, false, s12,
409  LATITUDE | LONGITUDE | AZIMUTH,
410  lat2, lon2, azi2, t, t, t, t, t);
411  }
412
413  /**
414  * See the documentation for Geodesic::Direct.
415  **********************************************************************/
416  Math::real Direct(real lat1, real lon1, real azi1, real s12,
417  real& lat2, real& lon2, real& azi2, real& m12)
418  const {
419  real t;
420  return GenDirect(lat1, lon1, azi1, false, s12,
421  LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH,
422  lat2, lon2, azi2, t, m12, t, t, t);
423  }
424
425  /**
426  * See the documentation for Geodesic::Direct.
427  **********************************************************************/
428  Math::real Direct(real lat1, real lon1, real azi1, real s12,
429  real& lat2, real& lon2, real& azi2,
430  real& M12, real& M21)
431  const {
432  real t;
433  return GenDirect(lat1, lon1, azi1, false, s12,
434  LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE,
435  lat2, lon2, azi2, t, t, M12, M21, t);
436  }
437
438  /**
439  * See the documentation for Geodesic::Direct.
440  **********************************************************************/
441  Math::real Direct(real lat1, real lon1, real azi1, real s12,
442  real& lat2, real& lon2, real& azi2,
443  real& m12, real& M12, real& M21)
444  const {
445  real t;
446  return GenDirect(lat1, lon1, azi1, false, s12,
447  LATITUDE | LONGITUDE | AZIMUTH |
448  REDUCEDLENGTH | GEODESICSCALE,
449  lat2, lon2, azi2, t, m12, M12, M21, t);
450  }
451  ///@}
452
453  /** \name Direct geodesic problem specified in terms of arc length.
454  **********************************************************************/
455  ///@{
456  /**
457  * Solve the direct geodesic problem where the length of the geodesic
458  * is specified in terms of arc length.
459  *
460  * @param[in] lat1 latitude of point 1 (degrees).
461  * @param[in] lon1 longitude of point 1 (degrees).
462  * @param[in] azi1 azimuth at point 1 (degrees).
463  * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
464  * be negative.
465  * @param[out] lat2 latitude of point 2 (degrees).
466  * @param[out] lon2 longitude of point 2 (degrees).
467  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
468  * @param[out] s12 distance between point 1 and point 2 (meters).
469  * @param[out] m12 reduced length of geodesic (meters).
470  * @param[out] M12 geodesic scale of point 2 relative to point 1
471  * (dimensionless).
472  * @param[out] M21 geodesic scale of point 1 relative to point 2
473  * (dimensionless).
474  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
475  *
476  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The values of
477  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
478  * 180&deg;).
479  *
480  * If either point is at a pole, the azimuth is defined by keeping the
481  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
482  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
483  * 180&deg; signifies a geodesic which is not a shortest path. (For a
484  * prolate ellipsoid, an additional condition is necessary for a shortest
485  * path: the longitudinal extent must not exceed of 180&deg;.)
486  *
487  * The following functions are overloaded versions of Geodesic::Direct
488  * which omit some of the output parameters.
489  **********************************************************************/
490  void ArcDirect(real lat1, real lon1, real azi1, real a12,
491  real& lat2, real& lon2, real& azi2, real& s12,
492  real& m12, real& M12, real& M21, real& S12)
493  const {
494  GenDirect(lat1, lon1, azi1, true, a12,
495  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
496  REDUCEDLENGTH | GEODESICSCALE | AREA,
497  lat2, lon2, azi2, s12, m12, M12, M21, S12);
498  }
499
500  /**
501  * See the documentation for Geodesic::ArcDirect.
502  **********************************************************************/
503  void ArcDirect(real lat1, real lon1, real azi1, real a12,
504  real& lat2, real& lon2) const {
505  real t;
506  GenDirect(lat1, lon1, azi1, true, a12,
507  LATITUDE | LONGITUDE,
508  lat2, lon2, t, t, t, t, t, t);
509  }
510
511  /**
512  * See the documentation for Geodesic::ArcDirect.
513  **********************************************************************/
514  void ArcDirect(real lat1, real lon1, real azi1, real a12,
515  real& lat2, real& lon2, real& azi2) const {
516  real t;
517  GenDirect(lat1, lon1, azi1, true, a12,
518  LATITUDE | LONGITUDE | AZIMUTH,
519  lat2, lon2, azi2, t, t, t, t, t);
520  }
521
522  /**
523  * See the documentation for Geodesic::ArcDirect.
524  **********************************************************************/
525  void ArcDirect(real lat1, real lon1, real azi1, real a12,
526  real& lat2, real& lon2, real& azi2, real& s12)
527  const {
528  real t;
529  GenDirect(lat1, lon1, azi1, true, a12,
530  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
531  lat2, lon2, azi2, s12, t, t, t, t);
532  }
533
534  /**
535  * See the documentation for Geodesic::ArcDirect.
536  **********************************************************************/
537  void ArcDirect(real lat1, real lon1, real azi1, real a12,
538  real& lat2, real& lon2, real& azi2,
539  real& s12, real& m12) const {
540  real t;
541  GenDirect(lat1, lon1, azi1, true, a12,
542  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
543  REDUCEDLENGTH,
544  lat2, lon2, azi2, s12, m12, t, t, t);
545  }
546
547  /**
548  * See the documentation for Geodesic::ArcDirect.
549  **********************************************************************/
550  void ArcDirect(real lat1, real lon1, real azi1, real a12,
551  real& lat2, real& lon2, real& azi2, real& s12,
552  real& M12, real& M21) const {
553  real t;
554  GenDirect(lat1, lon1, azi1, true, a12,
555  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
556  GEODESICSCALE,
557  lat2, lon2, azi2, s12, t, M12, M21, t);
558  }
559
560  /**
561  * See the documentation for Geodesic::ArcDirect.
562  **********************************************************************/
563  void ArcDirect(real lat1, real lon1, real azi1, real a12,
564  real& lat2, real& lon2, real& azi2, real& s12,
565  real& m12, real& M12, real& M21) const {
566  real t;
567  GenDirect(lat1, lon1, azi1, true, a12,
568  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
569  REDUCEDLENGTH | GEODESICSCALE,
570  lat2, lon2, azi2, s12, m12, M12, M21, t);
571  }
572  ///@}
573
574  /** \name General version of the direct geodesic solution.
575  **********************************************************************/
576  ///@{
577
578  /**
579  * The general direct geodesic problem. Geodesic::Direct and
580  * Geodesic::ArcDirect are defined in terms of this function.
581  *
582  * @param[in] lat1 latitude of point 1 (degrees).
583  * @param[in] lon1 longitude of point 1 (degrees).
584  * @param[in] azi1 azimuth at point 1 (degrees).
585  * @param[in] arcmode boolean flag determining the meaning of the \e
586  * s12_a12.
587  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
588  * point 1 and point 2 (meters); otherwise it is the arc length between
589  * point 1 and point 2 (degrees); it can be negative.
590  * @param[in] outmask a bitor'ed combination of Geodesic::mask values
591  * specifying which of the following parameters should be set.
592  * @param[out] lat2 latitude of point 2 (degrees).
593  * @param[out] lon2 longitude of point 2 (degrees).
594  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
595  * @param[out] s12 distance between point 1 and point 2 (meters).
596  * @param[out] m12 reduced length of geodesic (meters).
597  * @param[out] M12 geodesic scale of point 2 relative to point 1
598  * (dimensionless).
599  * @param[out] M21 geodesic scale of point 1 relative to point 2
600  * (dimensionless).
601  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
602  * @return \e a12 arc length of between point 1 and point 2 (degrees).
603  *
604  * The Geodesic::mask values possible for \e outmask are
605  * - \e outmask |= Geodesic::LATITUDE for the latitude \e lat2;
606  * - \e outmask |= Geodesic::LONGITUDE for the latitude \e lon2;
607  * - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2;
608  * - \e outmask |= Geodesic::DISTANCE for the distance \e s12;
609  * - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e
610  * m12;
611  * - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e
612  * M12 and \e M21;
613  * - \e outmask |= Geodesic::AREA for the area \e S12;
614  * - \e outmask |= Geodesic::ALL for all of the above;
615  * - \e outmask |= Geodesic::LONG_UNROLL to unroll \e lon2 instead of
616  * wrapping it into the range [&minus;180&deg;, 180&deg;).
617  * .
618  * The function value \e a12 is always computed and returned and this
619  * equals \e s12_a12 is \e arcmode is true. If \e outmask includes
620  * Geodesic::DISTANCE and \e arcmode is false, then \e s12 = \e s12_a12.
621  * It is not necessary to include Geodesic::DISTANCE_IN in \e outmask; this
622  * is automatically included is \e arcmode is false.
623  *
624  * With the Geodesic::LONG_UNROLL bit set, the quantity \e lon2 &minus; \e
625  * lon1 indicates how many times and in what sense the geodesic encircles
626  * the ellipsoid.
627  **********************************************************************/
628  Math::real GenDirect(real lat1, real lon1, real azi1,
629  bool arcmode, real s12_a12, unsigned outmask,
630  real& lat2, real& lon2, real& azi2,
631  real& s12, real& m12, real& M12, real& M21,
632  real& S12) const;
633  ///@}
634
635  /** \name Inverse geodesic problem.
636  **********************************************************************/
637  ///@{
638  /**
639  * Solve the inverse geodesic problem.
640  *
641  * @param[in] lat1 latitude of point 1 (degrees).
642  * @param[in] lon1 longitude of point 1 (degrees).
643  * @param[in] lat2 latitude of point 2 (degrees).
644  * @param[in] lon2 longitude of point 2 (degrees).
645  * @param[out] s12 distance between point 1 and point 2 (meters).
646  * @param[out] azi1 azimuth at point 1 (degrees).
647  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
648  * @param[out] m12 reduced length of geodesic (meters).
649  * @param[out] M12 geodesic scale of point 2 relative to point 1
650  * (dimensionless).
651  * @param[out] M21 geodesic scale of point 1 relative to point 2
652  * (dimensionless).
653  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
654  * @return \e a12 arc length of between point 1 and point 2 (degrees).
655  *
656  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
657  * The values of \e azi1 and \e azi2 returned are in the range
658  * [&minus;180&deg;, 180&deg;).
659  *
660  * If either point is at a pole, the azimuth is defined by keeping the
661  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
662  * and taking the limit &epsilon; &rarr; 0+.
663  *
664  * The solution to the inverse problem is found using Newton's method. If
665  * this fails to converge (this is very unlikely in geodetic applications
666  * but does occur for very eccentric ellipsoids), then the bisection method
667  * is used to refine the solution.
668  *
669  * The following functions are overloaded versions of Geodesic::Inverse
670  * which omit some of the output parameters. Note, however, that the arc
671  * length is always computed and returned as the function value.
672  **********************************************************************/
673  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
674  real& s12, real& azi1, real& azi2, real& m12,
675  real& M12, real& M21, real& S12) const {
676  return GenInverse(lat1, lon1, lat2, lon2,
677  DISTANCE | AZIMUTH |
678  REDUCEDLENGTH | GEODESICSCALE | AREA,
679  s12, azi1, azi2, m12, M12, M21, S12);
680  }
681
682  /**
683  * See the documentation for Geodesic::Inverse.
684  **********************************************************************/
685  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
686  real& s12) const {
687  real t;
688  return GenInverse(lat1, lon1, lat2, lon2,
689  DISTANCE,
690  s12, t, t, t, t, t, t);
691  }
692
693  /**
694  * See the documentation for Geodesic::Inverse.
695  **********************************************************************/
696  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
697  real& azi1, real& azi2) const {
698  real t;
699  return GenInverse(lat1, lon1, lat2, lon2,
700  AZIMUTH,
701  t, azi1, azi2, t, t, t, t);
702  }
703
704  /**
705  * See the documentation for Geodesic::Inverse.
706  **********************************************************************/
707  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
708  real& s12, real& azi1, real& azi2)
709  const {
710  real t;
711  return GenInverse(lat1, lon1, lat2, lon2,
712  DISTANCE | AZIMUTH,
713  s12, azi1, azi2, t, t, t, t);
714  }
715
716  /**
717  * See the documentation for Geodesic::Inverse.
718  **********************************************************************/
719  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
720  real& s12, real& azi1, real& azi2, real& m12)
721  const {
722  real t;
723  return GenInverse(lat1, lon1, lat2, lon2,
724  DISTANCE | AZIMUTH | REDUCEDLENGTH,
725  s12, azi1, azi2, m12, t, t, t);
726  }
727
728  /**
729  * See the documentation for Geodesic::Inverse.
730  **********************************************************************/
731  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
732  real& s12, real& azi1, real& azi2,
733  real& M12, real& M21) const {
734  real t;
735  return GenInverse(lat1, lon1, lat2, lon2,
736  DISTANCE | AZIMUTH | GEODESICSCALE,
737  s12, azi1, azi2, t, M12, M21, t);
738  }
739
740  /**
741  * See the documentation for Geodesic::Inverse.
742  **********************************************************************/
743  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
744  real& s12, real& azi1, real& azi2, real& m12,
745  real& M12, real& M21) const {
746  real t;
747  return GenInverse(lat1, lon1, lat2, lon2,
748  DISTANCE | AZIMUTH |
749  REDUCEDLENGTH | GEODESICSCALE,
750  s12, azi1, azi2, m12, M12, M21, t);
751  }
752  ///@}
753
754  /** \name General version of inverse geodesic solution.
755  **********************************************************************/
756  ///@{
757  /**
758  * The general inverse geodesic calculation. Geodesic::Inverse is defined
759  * in terms of this function.
760  *
761  * @param[in] lat1 latitude of point 1 (degrees).
762  * @param[in] lon1 longitude of point 1 (degrees).
763  * @param[in] lat2 latitude of point 2 (degrees).
764  * @param[in] lon2 longitude of point 2 (degrees).
765  * @param[in] outmask a bitor'ed combination of Geodesic::mask values
766  * specifying which of the following parameters should be set.
767  * @param[out] s12 distance between point 1 and point 2 (meters).
768  * @param[out] azi1 azimuth at point 1 (degrees).
769  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
770  * @param[out] m12 reduced length of geodesic (meters).
771  * @param[out] M12 geodesic scale of point 2 relative to point 1
772  * (dimensionless).
773  * @param[out] M21 geodesic scale of point 1 relative to point 2
774  * (dimensionless).
775  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
776  * @return \e a12 arc length of between point 1 and point 2 (degrees).
777  *
778  * The Geodesic::mask values possible for \e outmask are
779  * - \e outmask |= Geodesic::DISTANCE for the distance \e s12;
780  * - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2;
781  * - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e
782  * m12;
783  * - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e
784  * M12 and \e M21;
785  * - \e outmask |= Geodesic::AREA for the area \e S12;
786  * - \e outmask |= Geodesic::ALL for all of the above.
787  * .
788  * The arc length is always computed and returned as the function value.
789  **********************************************************************/
790  Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
792  real& s12, real& azi1, real& azi2,
793  real& m12, real& M12, real& M21, real& S12) const;
794  ///@}
795
796  /** \name Interface to GeodesicLine.
797  **********************************************************************/
798  ///@{
799
800  /**
801  * Set up to compute several points on a single geodesic.
802  *
803  * @param[in] lat1 latitude of point 1 (degrees).
804  * @param[in] lon1 longitude of point 1 (degrees).
805  * @param[in] azi1 azimuth at point 1 (degrees).
806  * @param[in] caps bitor'ed combination of Geodesic::mask values
807  * specifying the capabilities the GeodesicLine object should possess,
808  * i.e., which quantities can be returned in calls to
809  * GeodesicLine::Position.
810  * @return a GeodesicLine object.
811  *
812  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
813  *
814  * The Geodesic::mask values are
815  * - \e caps |= Geodesic::LATITUDE for the latitude \e lat2; this is
816  * added automatically;
817  * - \e caps |= Geodesic::LONGITUDE for the latitude \e lon2;
818  * - \e caps |= Geodesic::AZIMUTH for the azimuth \e azi2; this is
819  * added automatically;
820  * - \e caps |= Geodesic::DISTANCE for the distance \e s12;
821  * - \e caps |= Geodesic::REDUCEDLENGTH for the reduced length \e m12;
822  * - \e caps |= Geodesic::GEODESICSCALE for the geodesic scales \e M12
823  * and \e M21;
824  * - \e caps |= Geodesic::AREA for the area \e S12;
825  * - \e caps |= Geodesic::DISTANCE_IN permits the length of the
826  * geodesic to be given in terms of \e s12; without this capability the
827  * length can only be specified in terms of arc length;
828  * - \e caps |= Geodesic::ALL for all of the above.
829  * .
830  * The default value of \e caps is Geodesic::ALL.
831  *
832  * If the point is at a pole, the azimuth is defined by keeping \e lon1
833  * fixed, writing \e lat1 = &plusmn;(90 &minus; &epsilon;), and taking the
834  * limit &epsilon; &rarr; 0+.
835  **********************************************************************/
836  GeodesicLine Line(real lat1, real lon1, real azi1, unsigned caps = ALL)
837  const;
838
839  /**
840  * Define a GeodesicLine in terms of the inverse geodesic problem.
841  *
842  * @param[in] lat1 latitude of point 1 (degrees).
843  * @param[in] lon1 longitude of point 1 (degrees).
844  * @param[in] lat2 latitude of point 2 (degrees).
845  * @param[in] lon2 longitude of point 2 (degrees).
846  * @param[in] caps bitor'ed combination of Geodesic::mask values
847  * specifying the capabilities the GeodesicLine object should possess,
848  * i.e., which quantities can be returned in calls to
849  * GeodesicLine::Position.
850  * @return a GeodesicLine object.
851  *
852  * This function sets point 3 of the GeodesicLine to correspond to point 2
853  * of the inverse geodesic problem.
854  *
855  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
856  **********************************************************************/
857  GeodesicLine InverseLine(real lat1, real lon1, real lat2, real lon2,
858  unsigned caps = ALL) const;
859
860  /**
861  * Define a GeodesicLine in terms of the direct geodesic problem specified
862  * in terms of distance.
863  *
864  * @param[in] lat1 latitude of point 1 (degrees).
865  * @param[in] lon1 longitude of point 1 (degrees).
866  * @param[in] azi1 azimuth at point 1 (degrees).
867  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
868  * negative.
869  * @param[in] caps bitor'ed combination of Geodesic::mask values
870  * specifying the capabilities the GeodesicLine object should possess,
871  * i.e., which quantities can be returned in calls to
872  * GeodesicLine::Position.
873  * @return a GeodesicLine object.
874  *
875  * This function sets point 3 of the GeodesicLine to correspond to point 2
876  * of the direct geodesic problem.
877  *
878  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
879  **********************************************************************/
880  GeodesicLine DirectLine(real lat1, real lon1, real azi1, real s12,
881  unsigned caps = ALL) const;
882
883  /**
884  * Define a GeodesicLine in terms of the direct geodesic problem specified
885  * in terms of arc length.
886  *
887  * @param[in] lat1 latitude of point 1 (degrees).
888  * @param[in] lon1 longitude of point 1 (degrees).
889  * @param[in] azi1 azimuth at point 1 (degrees).
890  * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
891  * be negative.
892  * @param[in] caps bitor'ed combination of Geodesic::mask values
893  * specifying the capabilities the GeodesicLine object should possess,
894  * i.e., which quantities can be returned in calls to
895  * GeodesicLine::Position.
896  * @return a GeodesicLine object.
897  *
898  * This function sets point 3 of the GeodesicLine to correspond to point 2
899  * of the direct geodesic problem.
900  *
901  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
902  **********************************************************************/
903  GeodesicLine ArcDirectLine(real lat1, real lon1, real azi1, real a12,
904  unsigned caps = ALL) const;
905
906  /**
907  * Define a GeodesicLine in terms of the direct geodesic problem specified
908  * in terms of either distance or arc length.
909  *
910  * @param[in] lat1 latitude of point 1 (degrees).
911  * @param[in] lon1 longitude of point 1 (degrees).
912  * @param[in] azi1 azimuth at point 1 (degrees).
913  * @param[in] arcmode boolean flag determining the meaning of the \e
914  * s12_a12.
915  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
916  * point 1 and point 2 (meters); otherwise it is the arc length between
917  * point 1 and point 2 (degrees); it can be negative.
918  * @param[in] caps bitor'ed combination of Geodesic::mask values
919  * specifying the capabilities the GeodesicLine object should possess,
920  * i.e., which quantities can be returned in calls to
921  * GeodesicLine::Position.
922  * @return a GeodesicLine object.
923  *
924  * This function sets point 3 of the GeodesicLine to correspond to point 2
925  * of the direct geodesic problem.
926  *
927  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
928  **********************************************************************/
929  GeodesicLine GenDirectLine(real lat1, real lon1, real azi1,
930  bool arcmode, real s12_a12,
931  unsigned caps = ALL) const;
932  ///@}
933
934  /** \name Inspector functions.
935  **********************************************************************/
936  ///@{
937
938  /**
939  * @return \e a the equatorial radius of the ellipsoid (meters). This is
940  * the value used in the constructor.
941  **********************************************************************/
942  Math::real MajorRadius() const { return _a; }
943
944  /**
945  * @return \e f the flattening of the ellipsoid. This is the
946  * value used in the constructor.
947  **********************************************************************/
948  Math::real Flattening() const { return _f; }
949
950  /**
951  * @return total area of ellipsoid in meters<sup>2</sup>. The area of a
952  * polygon encircling a pole can be found by adding
953  * Geodesic::EllipsoidArea()/2 to the sum of \e S12 for each side of the
954  * polygon.
955  **********************************************************************/
957  { return 4 * Math::pi() * _c2; }
958  ///@}
959
960  /**
961  * A global instantiation of Geodesic with the parameters for the WGS84
962  * ellipsoid.
963  **********************************************************************/
964  static const Geodesic& WGS84();
965
966  };
967
968 } // namespace GeographicLib
969
970 #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:550
static T pi()
Definition: Math.hpp:202
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:90
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &m12) const
Definition: Geodesic.hpp:537
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:378
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12) const
Definition: Geodesic.hpp:416
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
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:673
Definition: Geodesic.hpp:942
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2) const
Definition: Geodesic.hpp:503
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:563
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2) const
Definition: Geodesic.hpp:392
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:731
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12) const
Definition: Geodesic.hpp:719
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:441
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2) const
Definition: Geodesic.hpp:707
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12) const
Definition: Geodesic.hpp:685
Math::real EllipsoidArea() const
Definition: Geodesic.hpp:956
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2) const
Definition: Geodesic.hpp:514
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:490
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:743
Header for GeographicLib::Constants class.
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &azi1, real &azi2) const
Definition: Geodesic.hpp:696
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12) const
Definition: Geodesic.hpp:525
#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:404
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:428
Math::real Flattening() const
Definition: Geodesic.hpp:948
Geodesic calculations
Definition: Geodesic.hpp:171