GeographicLib  2.0
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-2022) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
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; [&minus;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_UNROLL
101  };
102 
103  static real CosSeries(real sinx, real cosx, const real c[], int n);
104  static real Astroid(real x, real y);
105 
106  real _a, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2;
107  real _cC4x[nC4x_];
108 
109  void Lengths(const EllipticFunction& E,
110  real sig12,
111  real ssig1, real csig1, real dn1,
112  real ssig2, real csig2, real dn2,
113  real cbet1, real cbet2, unsigned outmask,
114  real& s12s, real& m12a, real& m0,
115  real& M12, real& M21) const;
116  real InverseStart(EllipticFunction& E,
117  real sbet1, real cbet1, real dn1,
118  real sbet2, real cbet2, real dn2,
119  real lam12, real slam12, real clam12,
120  real& salp1, real& calp1,
121  real& salp2, real& calp2, real& dnm) const;
122  real Lambda12(real sbet1, real cbet1, real dn1,
123  real sbet2, real cbet2, real dn2,
124  real salp1, real calp1, real slam120, real clam120,
125  real& salp2, real& calp2, real& sig12,
126  real& ssig1, real& csig1, real& ssig2, real& csig2,
127  EllipticFunction& E,
128  real& domg12, bool diffp, real& dlam12) const;
129  real GenInverse(real lat1, real lon1, real lat2, real lon2,
130  unsigned outmask, real& s12,
131  real& salp1, real& calp1, real& salp2, real& calp2,
132  real& m12, real& M12, real& M21, real& S12) const;
133 
134  // These are Maxima generated functions to provide series approximations to
135  // the integrals for the area.
136  void C4coeff();
137  void C4f(real k2, real c[]) const;
138  // Large coefficients are split so that lo contains the low 52 bits and hi
139  // the rest. This choice avoids double rounding with doubles and higher
140  // precision types. float coefficients will suffer double rounding;
141  // however the accuracy is already lousy for floats.
142  static Math::real reale(long long y, long long z) {
143  using std::ldexp;
144  return ldexp(real(y), 52) + z;
145  }
146 #if GEOGRAPHICLIB_GEODESICEXACT_ORDER > 30
147  // These are currently unused extended versions of reale needed really
148  // large coefficients (when using 64th order series). Such coefficients
149  // would overflow floats.
150  static Math::real reale(long long x, long long y, long long z) {
151  using std::ldexp;
152  return ldexp(real(x), 2*52) + (ldexp(real(y), 52) + z);
153  }
154  static Math::real reale(long long w, long long x,
155  long long y, long long z) {
156  using std::ldexp;
157  return ldexp(real(w), 3*52) +
158  (ldexp(real(x), 2*52) + (ldexp(real(y), 52) + z));
159  }
160 #endif
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  * A combination of the common capabilities: GeodesicExact::LATITUDE,
204  * GeodesicExact::LONGITUDE, GeodesicExact::AZIMUTH,
205  * GeodesicExact::DISTANCE.
206  * @hideinitializer
207  **********************************************************************/
208  STANDARD = LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
209  /**
210  * Allow distance \e s12 to be used as input in the direct geodesic
211  * problem.
212  * @hideinitializer
213  **********************************************************************/
214  DISTANCE_IN = 1U<<11 | CAP_E,
215  /**
216  * Calculate reduced length \e m12.
217  * @hideinitializer
218  **********************************************************************/
219  REDUCEDLENGTH = 1U<<12 | CAP_D,
220  /**
221  * Calculate geodesic scales \e M12 and \e M21.
222  * @hideinitializer
223  **********************************************************************/
224  GEODESICSCALE = 1U<<13 | CAP_D,
225  /**
226  * Calculate area \e S12.
227  * @hideinitializer
228  **********************************************************************/
229  AREA = 1U<<14 | CAP_C4,
230  /**
231  * Unroll \e lon2 in the direct calculation.
232  * @hideinitializer
233  **********************************************************************/
234  LONG_UNROLL = 1U<<15,
235  /**
236  * All capabilities, calculate everything. (GeodesicExact::LONG_UNROLL
237  * is not included in this mask.)
238  * @hideinitializer
239  **********************************************************************/
240  ALL = OUT_ALL| CAP_ALL,
241  };
242 
243  /** \name Constructor
244  **********************************************************************/
245  ///@{
246  /**
247  * Constructor for an ellipsoid with
248  *
249  * @param[in] a equatorial radius (meters).
250  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
251  * Negative \e f gives a prolate ellipsoid.
252  * @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
253  * positive.
254  **********************************************************************/
255  GeodesicExact(real a, real f);
256  ///@}
257 
258  /** \name Direct geodesic problem specified in terms of distance.
259  **********************************************************************/
260  ///@{
261  /**
262  * Perform the direct geodesic calculation where the length of the geodesic
263  * is specified in terms of distance.
264  *
265  * @param[in] lat1 latitude of point 1 (degrees).
266  * @param[in] lon1 longitude of point 1 (degrees).
267  * @param[in] azi1 azimuth at point 1 (degrees).
268  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
269  * signed.
270  * @param[out] lat2 latitude of point 2 (degrees).
271  * @param[out] lon2 longitude of point 2 (degrees).
272  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
273  * @param[out] m12 reduced length of geodesic (meters).
274  * @param[out] M12 geodesic scale of point 2 relative to point 1
275  * (dimensionless).
276  * @param[out] M21 geodesic scale of point 1 relative to point 2
277  * (dimensionless).
278  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
279  * @return \e a12 arc length of between point 1 and point 2 (degrees).
280  *
281  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The values of
282  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
283  * 180&deg;].
284  *
285  * If either point is at a pole, the azimuth is defined by keeping the
286  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
287  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
288  * 180&deg; signifies a geodesic which is not a shortest path. (For a
289  * prolate ellipsoid, an additional condition is necessary for a shortest
290  * path: the longitudinal extent must not exceed of 180&deg;.)
291  *
292  * The following functions are overloaded versions of GeodesicExact::Direct
293  * which omit some of the output parameters. Note, however, that the arc
294  * length is always computed and returned as the function value.
295  **********************************************************************/
296  Math::real Direct(real lat1, real lon1, real azi1, real s12,
297  real& lat2, real& lon2, real& azi2,
298  real& m12, real& M12, real& M21, real& S12)
299  const {
300  real t;
301  return GenDirect(lat1, lon1, azi1, false, s12,
302  LATITUDE | LONGITUDE | AZIMUTH |
303  REDUCEDLENGTH | GEODESICSCALE | AREA,
304  lat2, lon2, azi2, t, m12, M12, M21, S12);
305  }
306 
307  /**
308  * See the documentation for GeodesicExact::Direct.
309  **********************************************************************/
310  Math::real Direct(real lat1, real lon1, real azi1, real s12,
311  real& lat2, real& lon2)
312  const {
313  real t;
314  return GenDirect(lat1, lon1, azi1, false, s12,
315  LATITUDE | LONGITUDE,
316  lat2, lon2, t, t, t, t, t, t);
317  }
318 
319  /**
320  * See the documentation for GeodesicExact::Direct.
321  **********************************************************************/
322  Math::real Direct(real lat1, real lon1, real azi1, real s12,
323  real& lat2, real& lon2, real& azi2)
324  const {
325  real t;
326  return GenDirect(lat1, lon1, azi1, false, s12,
327  LATITUDE | LONGITUDE | AZIMUTH,
328  lat2, lon2, azi2, t, t, t, t, t);
329  }
330 
331  /**
332  * See the documentation for GeodesicExact::Direct.
333  **********************************************************************/
334  Math::real Direct(real lat1, real lon1, real azi1, real s12,
335  real& lat2, real& lon2, real& azi2, real& m12)
336  const {
337  real t;
338  return GenDirect(lat1, lon1, azi1, false, s12,
339  LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH,
340  lat2, lon2, azi2, t, m12, t, t, 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& M21)
349  const {
350  real t;
351  return GenDirect(lat1, lon1, azi1, false, s12,
352  LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE,
353  lat2, lon2, azi2, t, t, M12, M21, t);
354  }
355 
356  /**
357  * See the documentation for GeodesicExact::Direct.
358  **********************************************************************/
359  Math::real Direct(real lat1, real lon1, real azi1, real s12,
360  real& lat2, real& lon2, real& azi2,
361  real& m12, real& M12, real& M21)
362  const {
363  real t;
364  return GenDirect(lat1, lon1, azi1, false, s12,
365  LATITUDE | LONGITUDE | AZIMUTH |
366  REDUCEDLENGTH | GEODESICSCALE,
367  lat2, lon2, azi2, t, m12, M12, M21, t);
368  }
369  ///@}
370 
371  /** \name Direct geodesic problem specified in terms of arc length.
372  **********************************************************************/
373  ///@{
374  /**
375  * Perform the direct geodesic calculation where the length of the geodesic
376  * is specified in terms of arc length.
377  *
378  * @param[in] lat1 latitude of point 1 (degrees).
379  * @param[in] lon1 longitude of point 1 (degrees).
380  * @param[in] azi1 azimuth at point 1 (degrees).
381  * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
382  * be signed.
383  * @param[out] lat2 latitude of point 2 (degrees).
384  * @param[out] lon2 longitude of point 2 (degrees).
385  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
386  * @param[out] s12 distance between point 1 and point 2 (meters).
387  * @param[out] m12 reduced length of geodesic (meters).
388  * @param[out] M12 geodesic scale of point 2 relative to point 1
389  * (dimensionless).
390  * @param[out] M21 geodesic scale of point 1 relative to point 2
391  * (dimensionless).
392  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
393  *
394  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The values of
395  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
396  * 180&deg;].
397  *
398  * If either point is at a pole, the azimuth is defined by keeping the
399  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
400  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
401  * 180&deg; signifies a geodesic which is not a shortest path. (For a
402  * prolate ellipsoid, an additional condition is necessary for a shortest
403  * path: the longitudinal extent must not exceed of 180&deg;.)
404  *
405  * The following functions are overloaded versions of GeodesicExact::Direct
406  * which omit some of the output parameters.
407  **********************************************************************/
408  void ArcDirect(real lat1, real lon1, real azi1, real a12,
409  real& lat2, real& lon2, real& azi2, real& s12,
410  real& m12, real& M12, real& M21, real& S12)
411  const {
412  GenDirect(lat1, lon1, azi1, true, a12,
413  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
414  REDUCEDLENGTH | GEODESICSCALE | AREA,
415  lat2, lon2, azi2, s12, m12, M12, M21, S12);
416  }
417 
418  /**
419  * See the documentation for GeodesicExact::ArcDirect.
420  **********************************************************************/
421  void ArcDirect(real lat1, real lon1, real azi1, real a12,
422  real& lat2, real& lon2) const {
423  real t;
424  GenDirect(lat1, lon1, azi1, true, a12,
425  LATITUDE | LONGITUDE,
426  lat2, lon2, t, t, t, t, t, t);
427  }
428 
429  /**
430  * See the documentation for GeodesicExact::ArcDirect.
431  **********************************************************************/
432  void ArcDirect(real lat1, real lon1, real azi1, real a12,
433  real& lat2, real& lon2, real& azi2) const {
434  real t;
435  GenDirect(lat1, lon1, azi1, true, a12,
436  LATITUDE | LONGITUDE | AZIMUTH,
437  lat2, lon2, azi2, t, 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, real& s12)
445  const {
446  real t;
447  GenDirect(lat1, lon1, azi1, true, a12,
448  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
449  lat2, lon2, azi2, s12, t, 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,
457  real& s12, real& m12) const {
458  real t;
459  GenDirect(lat1, lon1, azi1, true, a12,
460  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
461  REDUCEDLENGTH,
462  lat2, lon2, azi2, s12, m12, t, t, 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& M21) const {
471  real t;
472  GenDirect(lat1, lon1, azi1, true, a12,
473  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
474  GEODESICSCALE,
475  lat2, lon2, azi2, s12, t, M12, M21, t);
476  }
477 
478  /**
479  * See the documentation for GeodesicExact::ArcDirect.
480  **********************************************************************/
481  void ArcDirect(real lat1, real lon1, real azi1, real a12,
482  real& lat2, real& lon2, real& azi2, real& s12,
483  real& m12, real& M12, real& M21) const {
484  real t;
485  GenDirect(lat1, lon1, azi1, true, a12,
486  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
487  REDUCEDLENGTH | GEODESICSCALE,
488  lat2, lon2, azi2, s12, m12, M12, M21, t);
489  }
490  ///@}
491 
492  /** \name General version of the direct geodesic solution.
493  **********************************************************************/
494  ///@{
495 
496  /**
497  * The general direct geodesic calculation. GeodesicExact::Direct and
498  * GeodesicExact::ArcDirect are defined in terms of this function.
499  *
500  * @param[in] lat1 latitude of point 1 (degrees).
501  * @param[in] lon1 longitude of point 1 (degrees).
502  * @param[in] azi1 azimuth at point 1 (degrees).
503  * @param[in] arcmode boolean flag determining the meaning of the second
504  * parameter.
505  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
506  * point 1 and point 2 (meters); otherwise it is the arc length between
507  * point 1 and point 2 (degrees); it can be signed.
508  * @param[in] outmask a bitor'ed combination of GeodesicExact::mask values
509  * specifying which of the following parameters should be set.
510  * @param[out] lat2 latitude of point 2 (degrees).
511  * @param[out] lon2 longitude of point 2 (degrees).
512  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
513  * @param[out] s12 distance between point 1 and point 2 (meters).
514  * @param[out] m12 reduced length of geodesic (meters).
515  * @param[out] M12 geodesic scale of point 2 relative to point 1
516  * (dimensionless).
517  * @param[out] M21 geodesic scale of point 1 relative to point 2
518  * (dimensionless).
519  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
520  * @return \e a12 arc length of between point 1 and point 2 (degrees).
521  *
522  * The GeodesicExact::mask values possible for \e outmask are
523  * - \e outmask |= GeodesicExact::LATITUDE for the latitude \e lat2;
524  * - \e outmask |= GeodesicExact::LONGITUDE for the latitude \e lon2;
525  * - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2;
526  * - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12;
527  * - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e
528  * m12;
529  * - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e
530  * M12 and \e M21;
531  * - \e outmask |= GeodesicExact::AREA for the area \e S12;
532  * - \e outmask |= GeodesicExact::ALL for all of the above;
533  * - \e outmask |= GeodesicExact::LONG_UNROLL to unroll \e lon2 instead of
534  * wrapping it into the range [&minus;180&deg;, 180&deg;].
535  * .
536  * The function value \e a12 is always computed and returned and this
537  * equals \e s12_a12 is \e arcmode is true. If \e outmask includes
538  * GeodesicExact::DISTANCE and \e arcmode is false, then \e s12 = \e
539  * s12_a12. It is not necessary to include GeodesicExact::DISTANCE_IN in
540  * \e outmask; this is automatically included is \e arcmode is false.
541  *
542  * With the GeodesicExact::LONG_UNROLL bit set, the quantity \e lon2
543  * &minus; \e lon1 indicates how many times and in what sense the geodesic
544  * encircles the ellipsoid.
545  **********************************************************************/
546  Math::real GenDirect(real lat1, real lon1, real azi1,
547  bool arcmode, real s12_a12, unsigned outmask,
548  real& lat2, real& lon2, real& azi2,
549  real& s12, real& m12, real& M12, real& M21,
550  real& S12) const;
551  ///@}
552 
553  /** \name Inverse geodesic problem.
554  **********************************************************************/
555  ///@{
556  /**
557  * Perform the inverse geodesic calculation.
558  *
559  * @param[in] lat1 latitude of point 1 (degrees).
560  * @param[in] lon1 longitude of point 1 (degrees).
561  * @param[in] lat2 latitude of point 2 (degrees).
562  * @param[in] lon2 longitude of point 2 (degrees).
563  * @param[out] s12 distance between point 1 and point 2 (meters).
564  * @param[out] azi1 azimuth at point 1 (degrees).
565  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
566  * @param[out] m12 reduced length of geodesic (meters).
567  * @param[out] M12 geodesic scale of point 2 relative to point 1
568  * (dimensionless).
569  * @param[out] M21 geodesic scale of point 1 relative to point 2
570  * (dimensionless).
571  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
572  * @return \e a12 arc length of between point 1 and point 2 (degrees).
573  *
574  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
575  * The values of \e azi1 and \e azi2 returned are in the range
576  * [&minus;180&deg;, 180&deg;].
577  *
578  * If either point is at a pole, the azimuth is defined by keeping the
579  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
580  * and taking the limit &epsilon; &rarr; 0+.
581  *
582  * The following functions are overloaded versions of
583  * GeodesicExact::Inverse which omit some of the output parameters. Note,
584  * however, that the arc length is always computed and returned as the
585  * function value.
586  **********************************************************************/
587  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
588  real& s12, real& azi1, real& azi2, real& m12,
589  real& M12, real& M21, real& S12) const {
590  return GenInverse(lat1, lon1, lat2, lon2,
591  DISTANCE | AZIMUTH |
592  REDUCEDLENGTH | GEODESICSCALE | AREA,
593  s12, azi1, azi2, m12, M12, M21, S12);
594  }
595 
596  /**
597  * See the documentation for GeodesicExact::Inverse.
598  **********************************************************************/
599  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
600  real& s12) const {
601  real t;
602  return GenInverse(lat1, lon1, lat2, lon2,
603  DISTANCE,
604  s12, t, t, t, t, t, t);
605  }
606 
607  /**
608  * See the documentation for GeodesicExact::Inverse.
609  **********************************************************************/
610  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
611  real& azi1, real& azi2) const {
612  real t;
613  return GenInverse(lat1, lon1, lat2, lon2,
614  AZIMUTH,
615  t, azi1, azi2, t, t, t, t);
616  }
617 
618  /**
619  * See the documentation for GeodesicExact::Inverse.
620  **********************************************************************/
621  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
622  real& s12, real& azi1, real& azi2)
623  const {
624  real t;
625  return GenInverse(lat1, lon1, lat2, lon2,
626  DISTANCE | AZIMUTH,
627  s12, azi1, azi2, t, t, t, t);
628  }
629 
630  /**
631  * See the documentation for GeodesicExact::Inverse.
632  **********************************************************************/
633  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
634  real& s12, real& azi1, real& azi2, real& m12)
635  const {
636  real t;
637  return GenInverse(lat1, lon1, lat2, lon2,
638  DISTANCE | AZIMUTH | REDUCEDLENGTH,
639  s12, azi1, azi2, m12, t, t, t);
640  }
641 
642  /**
643  * See the documentation for GeodesicExact::Inverse.
644  **********************************************************************/
645  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
646  real& s12, real& azi1, real& azi2,
647  real& M12, real& M21) const {
648  real t;
649  return GenInverse(lat1, lon1, lat2, lon2,
650  DISTANCE | AZIMUTH | GEODESICSCALE,
651  s12, azi1, azi2, t, M12, M21, t);
652  }
653 
654  /**
655  * See the documentation for GeodesicExact::Inverse.
656  **********************************************************************/
657  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
658  real& s12, real& azi1, real& azi2, real& m12,
659  real& M12, real& M21) const {
660  real t;
661  return GenInverse(lat1, lon1, lat2, lon2,
662  DISTANCE | AZIMUTH |
663  REDUCEDLENGTH | GEODESICSCALE,
664  s12, azi1, azi2, m12, M12, M21, t);
665  }
666  ///@}
667 
668  /** \name General version of inverse geodesic solution.
669  **********************************************************************/
670  ///@{
671  /**
672  * The general inverse geodesic calculation. GeodesicExact::Inverse is
673  * defined in terms of this function.
674  *
675  * @param[in] lat1 latitude of point 1 (degrees).
676  * @param[in] lon1 longitude of point 1 (degrees).
677  * @param[in] lat2 latitude of point 2 (degrees).
678  * @param[in] lon2 longitude of point 2 (degrees).
679  * @param[in] outmask a bitor'ed combination of GeodesicExact::mask values
680  * specifying which of the following parameters should be set.
681  * @param[out] s12 distance between point 1 and point 2 (meters).
682  * @param[out] azi1 azimuth at point 1 (degrees).
683  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
684  * @param[out] m12 reduced length of geodesic (meters).
685  * @param[out] M12 geodesic scale of point 2 relative to point 1
686  * (dimensionless).
687  * @param[out] M21 geodesic scale of point 1 relative to point 2
688  * (dimensionless).
689  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
690  * @return \e a12 arc length of between point 1 and point 2 (degrees).
691  *
692  * The GeodesicExact::mask values possible for \e outmask are
693  * - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12;
694  * - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2;
695  * - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e
696  * m12;
697  * - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e
698  * M12 and \e M21;
699  * - \e outmask |= GeodesicExact::AREA for the area \e S12;
700  * - \e outmask |= GeodesicExact::ALL for all of the above.
701  * .
702  * The arc length is always computed and returned as the function value.
703  **********************************************************************/
704  Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
705  unsigned outmask,
706  real& s12, real& azi1, real& azi2,
707  real& m12, real& M12, real& M21, real& S12) const;
708  ///@}
709 
710  /** \name Interface to GeodesicLineExact.
711  **********************************************************************/
712  ///@{
713 
714  /**
715  * Set up to compute several points on a single geodesic.
716  *
717  * @param[in] lat1 latitude of point 1 (degrees).
718  * @param[in] lon1 longitude of point 1 (degrees).
719  * @param[in] azi1 azimuth at point 1 (degrees).
720  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
721  * specifying the capabilities the GeodesicLineExact object should
722  * possess, i.e., which quantities can be returned in calls to
723  * GeodesicLineExact::Position.
724  * @return a GeodesicLineExact object.
725  *
726  * \e lat1 should be in the range [&minus;90&deg;, 90&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,
752  unsigned caps = ALL) const;
753 
754  /**
755  * Define a GeodesicLineExact in terms of the inverse geodesic problem.
756  *
757  * @param[in] lat1 latitude of point 1 (degrees).
758  * @param[in] lon1 longitude of point 1 (degrees).
759  * @param[in] lat2 latitude of point 2 (degrees).
760  * @param[in] lon2 longitude of point 2 (degrees).
761  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
762  * specifying the capabilities the GeodesicLineExact object should
763  * possess, i.e., which quantities can be returned in calls to
764  * GeodesicLineExact::Position.
765  * @return a GeodesicLineExact object.
766  *
767  * This function sets point 3 of the GeodesicLineExact to correspond to
768  * point 2 of the inverse geodesic problem.
769  *
770  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
771  **********************************************************************/
772  GeodesicLineExact InverseLine(real lat1, real lon1, real lat2, real lon2,
773  unsigned caps = ALL) const;
774 
775  /**
776  * Define a GeodesicLineExact in terms of the direct geodesic problem
777  * specified in terms of distance.
778  *
779  * @param[in] lat1 latitude of point 1 (degrees).
780  * @param[in] lon1 longitude of point 1 (degrees).
781  * @param[in] azi1 azimuth at point 1 (degrees).
782  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
783  * negative.
784  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
785  * specifying the capabilities the GeodesicLineExact object should
786  * possess, i.e., which quantities can be returned in calls to
787  * GeodesicLineExact::Position.
788  * @return a GeodesicLineExact object.
789  *
790  * This function sets point 3 of the GeodesicLineExact to correspond to
791  * point 2 of the direct geodesic problem.
792  *
793  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
794  **********************************************************************/
795  GeodesicLineExact DirectLine(real lat1, real lon1, real azi1, real s12,
796  unsigned caps = ALL) const;
797 
798  /**
799  * Define a GeodesicLineExact in terms of the direct geodesic problem
800  * specified in terms of arc length.
801  *
802  * @param[in] lat1 latitude of point 1 (degrees).
803  * @param[in] lon1 longitude of point 1 (degrees).
804  * @param[in] azi1 azimuth at point 1 (degrees).
805  * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
806  * be negative.
807  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
808  * specifying the capabilities the GeodesicLineExact object should
809  * possess, i.e., which quantities can be returned in calls to
810  * GeodesicLineExact::Position.
811  * @return a GeodesicLineExact object.
812  *
813  * This function sets point 3 of the GeodesicLineExact to correspond to
814  * point 2 of the direct geodesic problem.
815  *
816  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
817  **********************************************************************/
818  GeodesicLineExact ArcDirectLine(real lat1, real lon1, real azi1, real a12,
819  unsigned caps = ALL) const;
820 
821  /**
822  * Define a GeodesicLineExact in terms of the direct geodesic problem
823  * specified in terms of either distance or arc length.
824  *
825  * @param[in] lat1 latitude of point 1 (degrees).
826  * @param[in] lon1 longitude of point 1 (degrees).
827  * @param[in] azi1 azimuth at point 1 (degrees).
828  * @param[in] arcmode boolean flag determining the meaning of the \e
829  * s12_a12.
830  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
831  * point 1 and point 2 (meters); otherwise it is the arc length between
832  * point 1 and point 2 (degrees); it can be negative.
833  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
834  * specifying the capabilities the GeodesicLineExact object should
835  * possess, i.e., which quantities can be returned in calls to
836  * GeodesicLineExact::Position.
837  * @return a GeodesicLineExact object.
838  *
839  * This function sets point 3 of the GeodesicLineExact to correspond to
840  * point 2 of the direct geodesic problem.
841  *
842  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
843  **********************************************************************/
844  GeodesicLineExact GenDirectLine(real lat1, real lon1, real azi1,
845  bool arcmode, real s12_a12,
846  unsigned caps = ALL) const;
847  ///@}
848 
849  /** \name Inspector functions.
850  **********************************************************************/
851  ///@{
852 
853  /**
854  * @return \e a the equatorial radius of the ellipsoid (meters). This is
855  * the value used in the constructor.
856  **********************************************************************/
857  Math::real EquatorialRadius() const { return _a; }
858 
859  /**
860  * @return \e f the flattening of the ellipsoid. This is the
861  * value used in the constructor.
862  **********************************************************************/
863  Math::real Flattening() const { return _f; }
864 
865  /**
866  * @return total area of ellipsoid in meters<sup>2</sup>. The area of a
867  * polygon encircling a pole can be found by adding
868  * GeodesicExact::EllipsoidArea()/2 to the sum of \e S12 for each side of
869  * the polygon.
870  **********************************************************************/
872  { return 4 * Math::pi() * _c2; }
873  ///@}
874 
875  /**
876  * A global instantiation of GeodesicExact with the parameters for the
877  * WGS84 ellipsoid.
878  **********************************************************************/
879  static const GeodesicExact& WGS84();
880 
881  };
882 
883 } // namespace GeographicLib
884 
885 #endif // GEOGRAPHICLIB_GEODESICEXACT_HPP
Header for GeographicLib::Constants class.
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:67
Header for GeographicLib::EllipticFunction class.
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
#define GEOGRAPHICLIB_GEODESICEXACT_ORDER
Elliptic integrals and functions.
Exact geodesic calculations.
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, real &azi2, real &s12, real &M12, real &M21) const
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12) const
Math::real EllipsoidArea() const
Math::real Flattening() 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
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12) 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
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
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12) 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, real &m12, 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 &M21) const
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2) 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) const
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2) 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
Math::real EquatorialRadius() const
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2) const
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &m12) const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2) const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &M12, real &M21) const
static T pi()
Definition: Math.hpp:190
Namespace for GeographicLib.
Definition: Accumulator.cpp:12