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