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