GeographicLib  1.35
GeodesicExact.hpp
Go to the documentation of this file.
1 /**
2  * \file GeodesicExact.hpp
3  * \brief Header for GeographicLib::GeodesicExact class
4  *
5  * Copyright (c) Charles Karney (2012-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_GEODESICEXACT_HPP)
11 #define GEOGRAPHICLIB_GEODESICEXACT_HPP 1
12 
15 
16 #if !defined(GEOGRAPHICLIB_GEODESICEXACT_ORDER)
17 /**
18  * The order of the expansions used by GeodesicExact.
19  **********************************************************************/
20 # define GEOGRAPHICLIB_GEODESICEXACT_ORDER 30
21 #endif
22 
23 namespace GeographicLib {
24 
25  class GeodesicLineExact;
26 
27  /**
28  * \brief Exact geodesic calculations
29  *
30  * The equations for geodesics on an ellipsoid can be expressed in terms of
31  * incomplete elliptic integrals. The Geodesic class expands these integrals
32  * in a series in the flattening \e f and this provides an accurate solution
33  * for \e f &isin; [-0.01, 0.01]. The GeodesicExact class computes the
34  * ellitpic integrals directly and so provides a solution which is valid for
35  * all \e f. However, in practice, its use should be limited to about \e
36  * b/\e a &isin; [0.01, 100] or \e f &isin; [-99, 0.99].
37  *
38  * For the WGS84 ellipsoid, these classes are 2--3 times \e slower than the
39  * series solution and 2--3 times \e less \e accurate (because it's less easy
40  * to control round-off errors with the elliptic integral formulation); i.e.,
41  * the error is about 40 nm (40 nanometers) instead of 15 nm. However the
42  * error in the series solution scales as <i>f</i><sup>7</sup> while the
43  * error in the elliptic integral solution depends weakly on \e f. If the
44  * quarter meridian distance is 10000 km and the ratio \e b/\e a = 1 &minus;
45  * \e f is varied then the approximate maximum error (expressed as a
46  * distance) is <pre>
47  * 1 - f error (nm)
48  * 1/128 387
49  * 1/64 345
50  * 1/32 269
51  * 1/16 210
52  * 1/8 115
53  * 1/4 69
54  * 1/2 36
55  * 1 15
56  * 2 25
57  * 4 96
58  * 8 318
59  * 16 985
60  * 32 2352
61  * 64 6008
62  * 128 19024
63  * </pre>
64  *
65  * The computation of the area in these classes is via a 30th order series.
66  * This gives accurate results for \e b/\e a &isin; [1/2, 2]; the accuracy is
67  * about 8 decimal digits for \e b/\e a &isin; [1/4, 4].
68  *
69  * See \ref geodellip for the formulation. See the documentation on the
70  * Geodesic class for additional information on the geodesic problems.
71  *
72  * Example of use:
73  * \include example-GeodesicExact.cpp
74  *
75  * <a href="GeodSolve.1.html">GeodSolve</a> is a command-line utility
76  * providing access to the functionality of GeodesicExact and
77  * GeodesicLineExact (via the -E option).
78  **********************************************************************/
79 
81  private:
82  typedef Math::real real;
83  friend class GeodesicLineExact;
84  static const int nC4_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER;
85  static const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
86  static const unsigned maxit1_ = 20;
87  static const unsigned maxit2_ = maxit1_ +
88  std::numeric_limits<real>::digits + 10;
89 
90  static const real tiny_;
91  static const real tol0_;
92  static const real tol1_;
93  static const real tol2_;
94  static const real tolb_;
95  static const real xthresh_;
96 
97  enum captype {
98  CAP_NONE = 0U,
99  CAP_E = 1U<<0,
100  // Skip 1U<<1 for compatibility with Geodesic (not required)
101  CAP_D = 1U<<2,
102  CAP_H = 1U<<3,
103  CAP_C4 = 1U<<4,
104  CAP_ALL = 0x1FU,
105  OUT_ALL = 0x7F80U,
106  };
107 
108  static real CosSeries(real sinx, real cosx, const real c[], int n)
109  throw();
110  static inline real AngRound(real x) throw() {
111  // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
112  // for reals = 0.7 pm on the earth if x is an angle in degrees. (This
113  // is about 1000 times more resolution than we get with angles around 90
114  // degrees.) We use this to avoid having to deal with near singular
115  // cases when x is non-zero but tiny (e.g., 1.0e-200).
116  const real z = 1/real(16);
117  volatile real y = std::abs(x);
118  // The compiler mustn't "simplify" z - (z - y) to y
119  y = y < z ? z - (z - y) : y;
120  return x < 0 ? -y : y;
121  }
122  static inline void SinCosNorm(real& sinx, real& cosx) throw() {
123  real r = Math::hypot(sinx, cosx);
124  sinx /= r;
125  cosx /= r;
126  }
127  static real Astroid(real x, real y) throw();
128 
129  real _a, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2;
130  real _C4x[nC4x_];
131 
132  void Lengths(const EllipticFunction& E,
133  real sig12,
134  real ssig1, real csig1, real dn1,
135  real ssig2, real csig2, real dn2,
136  real cbet1, real cbet2,
137  real& s12s, real& m12a, real& m0,
138  bool scalep, real& M12, real& M21) const throw();
139  real InverseStart(EllipticFunction& E,
140  real sbet1, real cbet1, real dn1,
141  real sbet2, real cbet2, real dn2,
142  real lam12,
143  real& salp1, real& calp1,
144  real& salp2, real& calp2, real& dnm) const throw();
145  real Lambda12(real sbet1, real cbet1, real dn1,
146  real sbet2, real cbet2, real dn2,
147  real salp1, real calp1,
148  real& salp2, real& calp2, real& sig12,
149  real& ssig1, real& csig1, real& ssig2, real& csig2,
150  EllipticFunction& E,
151  real& omg12, bool diffp, real& dlam12)
152  const throw();
153 
154  // These are Maxima generated functions to provide series approximations to
155  // the integrals for the area.
156  void C4coeff() throw();
157  void C4f(real k2, real c[]) const throw();
158 
159  public:
160 
161  /**
162  * Bit masks for what calculations to do. These masks do double duty.
163  * They signify to the GeodesicLineExact::GeodesicLineExact constructor and
164  * to GeodesicExact::Line what capabilities should be included in the
165  * GeodesicLineExact object. They also specify which results to return in
166  * the general routines GeodesicExact::GenDirect and
167  * GeodesicExact::GenInverse routines. GeodesicLineExact::mask is a
168  * duplication of this enum.
169  **********************************************************************/
170  enum mask {
171  /**
172  * No capabilities, no output.
173  * @hideinitializer
174  **********************************************************************/
175  NONE = 0U,
176  /**
177  * Calculate latitude \e lat2. (It's not necessary to include this as a
178  * capability to GeodesicLineExact because this is included by default.)
179  * @hideinitializer
180  **********************************************************************/
181  LATITUDE = 1U<<7 | CAP_NONE,
182  /**
183  * Calculate longitude \e lon2.
184  * @hideinitializer
185  **********************************************************************/
186  LONGITUDE = 1U<<8 | CAP_H,
187  /**
188  * Calculate azimuths \e azi1 and \e azi2. (It's not necessary to
189  * include this as a capability to GeodesicLineExact because this is
190  * included by default.)
191  * @hideinitializer
192  **********************************************************************/
193  AZIMUTH = 1U<<9 | CAP_NONE,
194  /**
195  * Calculate distance \e s12.
196  * @hideinitializer
197  **********************************************************************/
198  DISTANCE = 1U<<10 | CAP_E,
199  /**
200  * Allow distance \e s12 to be used as input in the direct geodesic
201  * problem.
202  * @hideinitializer
203  **********************************************************************/
204  DISTANCE_IN = 1U<<11 | CAP_E,
205  /**
206  * Calculate reduced length \e m12.
207  * @hideinitializer
208  **********************************************************************/
209  REDUCEDLENGTH = 1U<<12 | CAP_D,
210  /**
211  * Calculate geodesic scales \e M12 and \e M21.
212  * @hideinitializer
213  **********************************************************************/
214  GEODESICSCALE = 1U<<13 | CAP_D,
215  /**
216  * Calculate area \e S12.
217  * @hideinitializer
218  **********************************************************************/
219  AREA = 1U<<14 | CAP_C4,
220  /**
221  * All capabilities, calculate everything.
222  * @hideinitializer
223  **********************************************************************/
224  ALL = OUT_ALL| CAP_ALL,
225  };
226 
227  /** \name Constructor
228  **********************************************************************/
229  ///@{
230  /**
231  * Constructor for a ellipsoid with
232  *
233  * @param[in] a equatorial radius (meters).
234  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
235  * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flattening
236  * to 1/\e f.
237  * @exception GeographicErr if \e a or (1 &minus; \e f ) \e a is not
238  * positive.
239  **********************************************************************/
240  GeodesicExact(real a, real f);
241  ///@}
242 
243  /** \name Direct geodesic problem specified in terms of distance.
244  **********************************************************************/
245  ///@{
246  /**
247  * Perform the direct geodesic calculation where the length of the geodesic
248  * is specified in terms of distance.
249  *
250  * @param[in] lat1 latitude of point 1 (degrees).
251  * @param[in] lon1 longitude of point 1 (degrees).
252  * @param[in] azi1 azimuth at point 1 (degrees).
253  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
254  * signed.
255  * @param[out] lat2 latitude of point 2 (degrees).
256  * @param[out] lon2 longitude of point 2 (degrees).
257  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
258  * @param[out] m12 reduced length of geodesic (meters).
259  * @param[out] M12 geodesic scale of point 2 relative to point 1
260  * (dimensionless).
261  * @param[out] M21 geodesic scale of point 1 relative to point 2
262  * (dimensionless).
263  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
264  * @return \e a12 arc length of between point 1 and point 2 (degrees).
265  *
266  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
267  * azi1 should be in the range [&minus;540&deg;, 540&deg;). The values of
268  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
269  * 180&deg;).
270  *
271  * If either point is at a pole, the azimuth is defined by keeping the
272  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
273  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
274  * 180&deg; signifies a geodesic which is not a shortest path. (For a
275  * prolate ellipsoid, an additional condition is necessary for a shortest
276  * path: the longitudinal extent must not exceed of 180&deg;.)
277  *
278  * The following functions are overloaded versions of GeodesicExact::Direct
279  * which omit some of the output parameters. Note, however, that the arc
280  * length is always computed and returned as the function value.
281  **********************************************************************/
282  Math::real Direct(real lat1, real lon1, real azi1, real s12,
283  real& lat2, real& lon2, real& azi2,
284  real& m12, real& M12, real& M21, real& S12)
285  const throw() {
286  real t;
287  return GenDirect(lat1, lon1, azi1, false, s12,
288  LATITUDE | LONGITUDE | AZIMUTH |
289  REDUCEDLENGTH | GEODESICSCALE | AREA,
290  lat2, lon2, azi2, t, m12, M12, M21, S12);
291  }
292 
293  /**
294  * See the documentation for GeodesicExact::Direct.
295  **********************************************************************/
296  Math::real Direct(real lat1, real lon1, real azi1, real s12,
297  real& lat2, real& lon2)
298  const throw() {
299  real t;
300  return GenDirect(lat1, lon1, azi1, false, s12,
301  LATITUDE | LONGITUDE,
302  lat2, lon2, t, t, t, t, t, t);
303  }
304 
305  /**
306  * See the documentation for GeodesicExact::Direct.
307  **********************************************************************/
308  Math::real Direct(real lat1, real lon1, real azi1, real s12,
309  real& lat2, real& lon2, real& azi2)
310  const throw() {
311  real t;
312  return GenDirect(lat1, lon1, azi1, false, s12,
313  LATITUDE | LONGITUDE | AZIMUTH,
314  lat2, lon2, azi2, t, t, t, t, t);
315  }
316 
317  /**
318  * See the documentation for GeodesicExact::Direct.
319  **********************************************************************/
320  Math::real Direct(real lat1, real lon1, real azi1, real s12,
321  real& lat2, real& lon2, real& azi2, real& m12)
322  const throw() {
323  real t;
324  return GenDirect(lat1, lon1, azi1, false, s12,
325  LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH,
326  lat2, lon2, azi2, t, m12, t, t, t);
327  }
328 
329  /**
330  * See the documentation for GeodesicExact::Direct.
331  **********************************************************************/
332  Math::real Direct(real lat1, real lon1, real azi1, real s12,
333  real& lat2, real& lon2, real& azi2,
334  real& M12, real& M21)
335  const throw() {
336  real t;
337  return GenDirect(lat1, lon1, azi1, false, s12,
338  LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE,
339  lat2, lon2, azi2, t, t, M12, M21, t);
340  }
341 
342  /**
343  * See the documentation for GeodesicExact::Direct.
344  **********************************************************************/
345  Math::real Direct(real lat1, real lon1, real azi1, real s12,
346  real& lat2, real& lon2, real& azi2,
347  real& m12, real& M12, real& M21)
348  const throw() {
349  real t;
350  return GenDirect(lat1, lon1, azi1, false, s12,
351  LATITUDE | LONGITUDE | AZIMUTH |
352  REDUCEDLENGTH | GEODESICSCALE,
353  lat2, lon2, azi2, t, m12, M12, M21, t);
354  }
355  ///@}
356 
357  /** \name Direct geodesic problem specified in terms of arc length.
358  **********************************************************************/
359  ///@{
360  /**
361  * Perform the direct geodesic calculation where the length of the geodesic
362  * is specified in terms of arc length.
363  *
364  * @param[in] lat1 latitude of point 1 (degrees).
365  * @param[in] lon1 longitude of point 1 (degrees).
366  * @param[in] azi1 azimuth at point 1 (degrees).
367  * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
368  * be signed.
369  * @param[out] lat2 latitude of point 2 (degrees).
370  * @param[out] lon2 longitude of point 2 (degrees).
371  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
372  * @param[out] s12 distance between point 1 and point 2 (meters).
373  * @param[out] m12 reduced length of geodesic (meters).
374  * @param[out] M12 geodesic scale of point 2 relative to point 1
375  * (dimensionless).
376  * @param[out] M21 geodesic scale of point 1 relative to point 2
377  * (dimensionless).
378  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
379  *
380  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
381  * azi1 should be in the range [&minus;540&deg;, 540&deg;). The values of
382  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
383  * 180&deg;).
384  *
385  * If either point is at a pole, the azimuth is defined by keeping the
386  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
387  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
388  * 180&deg; signifies a geodesic which is not a shortest path. (For a
389  * prolate ellipsoid, an additional condition is necessary for a shortest
390  * path: the longitudinal extent must not exceed of 180&deg;.)
391  *
392  * The following functions are overloaded versions of GeodesicExact::Direct
393  * which omit some of the output parameters.
394  **********************************************************************/
395  void ArcDirect(real lat1, real lon1, real azi1, real a12,
396  real& lat2, real& lon2, real& azi2, real& s12,
397  real& m12, real& M12, real& M21, real& S12)
398  const throw() {
399  GenDirect(lat1, lon1, azi1, true, a12,
400  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
401  REDUCEDLENGTH | GEODESICSCALE | AREA,
402  lat2, lon2, azi2, s12, m12, M12, M21, S12);
403  }
404 
405  /**
406  * See the documentation for GeodesicExact::ArcDirect.
407  **********************************************************************/
408  void ArcDirect(real lat1, real lon1, real azi1, real a12,
409  real& lat2, real& lon2) const throw() {
410  real t;
411  GenDirect(lat1, lon1, azi1, true, a12,
412  LATITUDE | LONGITUDE,
413  lat2, lon2, t, t, t, t, t, t);
414  }
415 
416  /**
417  * See the documentation for GeodesicExact::ArcDirect.
418  **********************************************************************/
419  void ArcDirect(real lat1, real lon1, real azi1, real a12,
420  real& lat2, real& lon2, real& azi2) const throw() {
421  real t;
422  GenDirect(lat1, lon1, azi1, true, a12,
423  LATITUDE | LONGITUDE | AZIMUTH,
424  lat2, lon2, azi2, t, t, t, t, t);
425  }
426 
427  /**
428  * See the documentation for GeodesicExact::ArcDirect.
429  **********************************************************************/
430  void ArcDirect(real lat1, real lon1, real azi1, real a12,
431  real& lat2, real& lon2, real& azi2, real& s12)
432  const throw() {
433  real t;
434  GenDirect(lat1, lon1, azi1, true, a12,
435  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
436  lat2, lon2, azi2, s12, t, t, t, t);
437  }
438 
439  /**
440  * See the documentation for GeodesicExact::ArcDirect.
441  **********************************************************************/
442  void ArcDirect(real lat1, real lon1, real azi1, real a12,
443  real& lat2, real& lon2, real& azi2,
444  real& s12, real& m12) const throw() {
445  real t;
446  GenDirect(lat1, lon1, azi1, true, a12,
447  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
448  REDUCEDLENGTH,
449  lat2, lon2, azi2, s12, m12, t, t, t);
450  }
451 
452  /**
453  * See the documentation for GeodesicExact::ArcDirect.
454  **********************************************************************/
455  void ArcDirect(real lat1, real lon1, real azi1, real a12,
456  real& lat2, real& lon2, real& azi2, real& s12,
457  real& M12, real& M21) const throw() {
458  real t;
459  GenDirect(lat1, lon1, azi1, true, a12,
460  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
461  GEODESICSCALE,
462  lat2, lon2, azi2, s12, t, M12, M21, t);
463  }
464 
465  /**
466  * See the documentation for GeodesicExact::ArcDirect.
467  **********************************************************************/
468  void ArcDirect(real lat1, real lon1, real azi1, real a12,
469  real& lat2, real& lon2, real& azi2, real& s12,
470  real& m12, real& M12, real& M21) const throw() {
471  real t;
472  GenDirect(lat1, lon1, azi1, true, a12,
473  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
474  REDUCEDLENGTH | GEODESICSCALE,
475  lat2, lon2, azi2, s12, m12, M12, M21, t);
476  }
477  ///@}
478 
479  /** \name General version of the direct geodesic solution.
480  **********************************************************************/
481  ///@{
482 
483  /**
484  * The general direct geodesic calculation. GeodesicExact::Direct and
485  * GeodesicExact::ArcDirect are defined in terms of this function.
486  *
487  * @param[in] lat1 latitude of point 1 (degrees).
488  * @param[in] lon1 longitude of point 1 (degrees).
489  * @param[in] azi1 azimuth at point 1 (degrees).
490  * @param[in] arcmode boolean flag determining the meaning of the second
491  * parameter.
492  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
493  * point 1 and point 2 (meters); otherwise it is the arc length between
494  * point 1 and point 2 (degrees); it can be signed.
495  * @param[in] outmask a bitor'ed combination of GeodesicExact::mask values
496  * specifying which of the following parameters should be set.
497  * @param[out] lat2 latitude of point 2 (degrees).
498  * @param[out] lon2 longitude of point 2 (degrees).
499  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
500  * @param[out] s12 distance between point 1 and point 2 (meters).
501  * @param[out] m12 reduced length of geodesic (meters).
502  * @param[out] M12 geodesic scale of point 2 relative to point 1
503  * (dimensionless).
504  * @param[out] M21 geodesic scale of point 1 relative to point 2
505  * (dimensionless).
506  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
507  * @return \e a12 arc length of between point 1 and point 2 (degrees).
508  *
509  * The GeodesicExact::mask values possible for \e outmask are
510  * - \e outmask |= GeodesicExact::LATITUDE for the latitude \e lat2;
511  * - \e outmask |= GeodesicExact::LONGITUDE for the latitude \e lon2;
512  * - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2;
513  * - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12;
514  * - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e
515  * m12;
516  * - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e
517  * M12 and \e M21;
518  * - \e outmask |= GeodesicExact::AREA for the area \e S12;
519  * - \e outmask |= GeodesicExact::ALL for all of the above.
520  * .
521  * The function value \e a12 is always computed and returned and this
522  * equals \e s12_a12 is \e arcmode is true. If \e outmask includes
523  * GeodesicExact::DISTANCE and \e arcmode is false, then \e s12 = \e
524  * s12_a12. It is not necessary to include GeodesicExact::DISTANCE_IN in
525  * \e outmask; this is automatically included is \e arcmode is false.
526  **********************************************************************/
527  Math::real GenDirect(real lat1, real lon1, real azi1,
528  bool arcmode, real s12_a12, unsigned outmask,
529  real& lat2, real& lon2, real& azi2,
530  real& s12, real& m12, real& M12, real& M21,
531  real& S12) const throw();
532  ///@}
533 
534  /** \name Inverse geodesic problem.
535  **********************************************************************/
536  ///@{
537  /**
538  * Perform the inverse geodesic calculation.
539  *
540  * @param[in] lat1 latitude of point 1 (degrees).
541  * @param[in] lon1 longitude of point 1 (degrees).
542  * @param[in] lat2 latitude of point 2 (degrees).
543  * @param[in] lon2 longitude of point 2 (degrees).
544  * @param[out] s12 distance between point 1 and point 2 (meters).
545  * @param[out] azi1 azimuth at point 1 (degrees).
546  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
547  * @param[out] m12 reduced length of geodesic (meters).
548  * @param[out] M12 geodesic scale of point 2 relative to point 1
549  * (dimensionless).
550  * @param[out] M21 geodesic scale of point 1 relative to point 2
551  * (dimensionless).
552  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
553  * @return \e a12 arc length of between point 1 and point 2 (degrees).
554  *
555  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e
556  * lon1 and \e lon2 should be in the range [&minus;540&deg;, 540&deg;).
557  * The values of \e azi1 and \e azi2 returned are in the range
558  * [&minus;180&deg;, 180&deg;).
559  *
560  * If either point is at a pole, the azimuth is defined by keeping the
561  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
562  * and taking the limit &epsilon; &rarr; 0+.
563  *
564  * The following functions are overloaded versions of GeodesicExact::Inverse
565  * which omit some of the output parameters. Note, however, that the arc
566  * length is always computed and returned as the function value.
567  **********************************************************************/
568  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
569  real& s12, real& azi1, real& azi2, real& m12,
570  real& M12, real& M21, real& S12) const throw() {
571  return GenInverse(lat1, lon1, lat2, lon2,
572  DISTANCE | AZIMUTH |
573  REDUCEDLENGTH | GEODESICSCALE | AREA,
574  s12, azi1, azi2, m12, M12, M21, S12);
575  }
576 
577  /**
578  * See the documentation for GeodesicExact::Inverse.
579  **********************************************************************/
580  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
581  real& s12) const throw() {
582  real t;
583  return GenInverse(lat1, lon1, lat2, lon2,
584  DISTANCE,
585  s12, t, t, t, t, t, t);
586  }
587 
588  /**
589  * See the documentation for GeodesicExact::Inverse.
590  **********************************************************************/
591  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
592  real& azi1, real& azi2) const throw() {
593  real t;
594  return GenInverse(lat1, lon1, lat2, lon2,
595  AZIMUTH,
596  t, azi1, azi2, t, t, t, t);
597  }
598 
599  /**
600  * See the documentation for GeodesicExact::Inverse.
601  **********************************************************************/
602  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
603  real& s12, real& azi1, real& azi2)
604  const throw() {
605  real t;
606  return GenInverse(lat1, lon1, lat2, lon2,
607  DISTANCE | AZIMUTH,
608  s12, azi1, azi2, t, t, t, t);
609  }
610 
611  /**
612  * See the documentation for GeodesicExact::Inverse.
613  **********************************************************************/
614  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
615  real& s12, real& azi1, real& azi2, real& m12)
616  const throw() {
617  real t;
618  return GenInverse(lat1, lon1, lat2, lon2,
619  DISTANCE | AZIMUTH | REDUCEDLENGTH,
620  s12, azi1, azi2, m12, t, t, t);
621  }
622 
623  /**
624  * See the documentation for GeodesicExact::Inverse.
625  **********************************************************************/
626  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
627  real& s12, real& azi1, real& azi2,
628  real& M12, real& M21) const throw() {
629  real t;
630  return GenInverse(lat1, lon1, lat2, lon2,
631  DISTANCE | AZIMUTH | GEODESICSCALE,
632  s12, azi1, azi2, t, M12, M21, t);
633  }
634 
635  /**
636  * See the documentation for GeodesicExact::Inverse.
637  **********************************************************************/
638  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
639  real& s12, real& azi1, real& azi2, real& m12,
640  real& M12, real& M21) const throw() {
641  real t;
642  return GenInverse(lat1, lon1, lat2, lon2,
643  DISTANCE | AZIMUTH |
644  REDUCEDLENGTH | GEODESICSCALE,
645  s12, azi1, azi2, m12, M12, M21, t);
646  }
647  ///@}
648 
649  /** \name General version of inverse geodesic solution.
650  **********************************************************************/
651  ///@{
652  /**
653  * The general inverse geodesic calculation. GeodesicExact::Inverse is
654  * defined in terms of this function.
655  *
656  * @param[in] lat1 latitude of point 1 (degrees).
657  * @param[in] lon1 longitude of point 1 (degrees).
658  * @param[in] lat2 latitude of point 2 (degrees).
659  * @param[in] lon2 longitude of point 2 (degrees).
660  * @param[in] outmask a bitor'ed combination of GeodesicExact::mask values
661  * specifying which of the following parameters should be set.
662  * @param[out] s12 distance between point 1 and point 2 (meters).
663  * @param[out] azi1 azimuth at point 1 (degrees).
664  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
665  * @param[out] m12 reduced length of geodesic (meters).
666  * @param[out] M12 geodesic scale of point 2 relative to point 1
667  * (dimensionless).
668  * @param[out] M21 geodesic scale of point 1 relative to point 2
669  * (dimensionless).
670  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
671  * @return \e a12 arc length of between point 1 and point 2 (degrees).
672  *
673  * The GeodesicExact::mask values possible for \e outmask are
674  * - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12;
675  * - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2;
676  * - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e
677  * m12;
678  * - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e
679  * M12 and \e M21;
680  * - \e outmask |= GeodesicExact::AREA for the area \e S12;
681  * - \e outmask |= GeodesicExact::ALL for all of the above.
682  * .
683  * The arc length is always computed and returned as the function value.
684  **********************************************************************/
685  Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
686  unsigned outmask,
687  real& s12, real& azi1, real& azi2,
688  real& m12, real& M12, real& M21, real& S12)
689  const throw();
690  ///@}
691 
692  /** \name Interface to GeodesicLineExact.
693  **********************************************************************/
694  ///@{
695 
696  /**
697  * Set up to compute several points on a single geodesic.
698  *
699  * @param[in] lat1 latitude of point 1 (degrees).
700  * @param[in] lon1 longitude of point 1 (degrees).
701  * @param[in] azi1 azimuth at point 1 (degrees).
702  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
703  * specifying the capabilities the GeodesicLineExact object should
704  * possess, i.e., which quantities can be returned in calls to
705  * GeodesicLineExact::Position.
706  * @return a GeodesicLineExact object.
707  *
708  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
709  * azi1 should be in the range [&minus;540&deg;, 540&deg;).
710  *
711  * The GeodesicExact::mask values are
712  * - \e caps |= GeodesicExact::LATITUDE for the latitude \e lat2; this is
713  * added automatically;
714  * - \e caps |= GeodesicExact::LONGITUDE for the latitude \e lon2;
715  * - \e caps |= GeodesicExact::AZIMUTH for the azimuth \e azi2; this is
716  * added automatically;
717  * - \e caps |= GeodesicExact::DISTANCE for the distance \e s12;
718  * - \e caps |= GeodesicExact::REDUCEDLENGTH for the reduced length \e m12;
719  * - \e caps |= GeodesicExact::GEODESICSCALE for the geodesic scales \e M12
720  * and \e M21;
721  * - \e caps |= GeodesicExact::AREA for the area \e S12;
722  * - \e caps |= GeodesicExact::DISTANCE_IN permits the length of the
723  * geodesic to be given in terms of \e s12; without this capability the
724  * length can only be specified in terms of arc length;
725  * - \e caps |= GeodesicExact::ALL for all of the above.
726  * .
727  * The default value of \e caps is GeodesicExact::ALL which turns on all
728  * the capabilities.
729  *
730  * If the point is at a pole, the azimuth is defined by keeping \e lon1
731  * fixed, writing \e lat1 = &plusmn;(90 &minus; &epsilon;), and taking the
732  * limit &epsilon; &rarr; 0+.
733  **********************************************************************/
734  GeodesicLineExact Line(real lat1, real lon1, real azi1, unsigned caps = ALL)
735  const throw();
736 
737  ///@}
738 
739  /** \name Inspector functions.
740  **********************************************************************/
741  ///@{
742 
743  /**
744  * @return \e a the equatorial radius of the ellipsoid (meters). This is
745  * the value used in the constructor.
746  **********************************************************************/
747  Math::real MajorRadius() const throw() { return _a; }
748 
749  /**
750  * @return \e f the flattening of the ellipsoid. This is the
751  * value used in the constructor.
752  **********************************************************************/
753  Math::real Flattening() const throw() { return _f; }
754 
755  /// \cond SKIP
756  /**
757  * <b>DEPRECATED</b>
758  * @return \e r the inverse flattening of the ellipsoid.
759  **********************************************************************/
760  Math::real InverseFlattening() const throw() { return 1/_f; }
761  /// \endcond
762 
763  /**
764  * @return total area of ellipsoid in meters<sup>2</sup>. The area of a
765  * polygon encircling a pole can be found by adding
766  * GeodesicExact::EllipsoidArea()/2 to the sum of \e S12 for each side of
767  * the polygon.
768  **********************************************************************/
769  Math::real EllipsoidArea() const throw()
770  { return 4 * Math::pi<real>() * _c2; }
771  ///@}
772 
773  /**
774  * A global instantiation of GeodesicExact with the parameters for the WGS84
775  * ellipsoid.
776  **********************************************************************/
777  static const GeodesicExact WGS84;
778 
779  };
780 
781 } // namespace GeographicLib
782 
783 #endif // GEOGRAPHICLIB_GEODESICEXACT_HPP
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
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:52
Math::real Flattening() const
Math::real EllipsoidArea() const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12, real &M12, real &M21) const
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2) const
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:73
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12, real &M12, real &M21) const
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2) const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &M12, real &M21) const
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
static const GeodesicExact WGS84
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12) const
#define GEOGRAPHICLIB_GEODESICEXACT_ORDER
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2) const
static T hypot(T x, T y)
Definition: Math.hpp:165
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &M12, real &M21) const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &azi1, real &azi2) const
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &m12) const
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &M12, real &M21) const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12) const
Header for GeographicLib::EllipticFunction class.
Exact geodesic calculations.
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2) const
Header for GeographicLib::Constants class.
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12) const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12) const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2) const
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
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