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