GeographicLib  2.1.2
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 #include <GeographicLib/DST.hpp>
16 #include <vector>
17 
18 namespace GeographicLib {
19 
20  class GEOGRAPHICLIB_EXPORT DST;
21  class GeodesicLineExact;
22 
23  /**
24  * \brief Exact geodesic calculations
25  *
26  * The equations for geodesics on an ellipsoid can be expressed in terms of
27  * incomplete elliptic integrals. The Geodesic class expands these integrals
28  * in a series in the flattening \e f and this provides an accurate solution
29  * for \e f &isin; [-0.01, 0.01]. The GeodesicExact class computes the
30  * ellitpic integrals directly and so provides a solution which is valid for
31  * all \e f. However, in practice, its use should be limited to about
32  * <i>b</i>/\e a &isin; [0.01, 100] or \e f &isin; [&minus;99, 0.99].
33  *
34  * For the WGS84 ellipsoid, these classes are 2--3 times \e slower than the
35  * series solution and 2--3 times \e less \e accurate (because it's less easy
36  * to control round-off errors with the elliptic integral formulation); i.e.,
37  * the error is about 40 nm (40 nanometers) instead of 15 nm. However the
38  * error in the series solution scales as <i>f</i><sup>7</sup> while the
39  * error in the elliptic integral solution depends weakly on \e f. If the
40  * quarter meridian distance is 10000 km and the ratio <i>b</i>/\e a = 1
41  * &minus; \e f is varied then the approximate maximum error (expressed as a
42  * distance) is <pre>
43  * 1 - f error (nm)
44  * 1/128 387
45  * 1/64 345
46  * 1/32 269
47  * 1/16 210
48  * 1/8 115
49  * 1/4 69
50  * 1/2 36
51  * 1 15
52  * 2 25
53  * 4 96
54  * 8 318
55  * 16 985
56  * 32 2352
57  * 64 6008
58  * 128 19024
59  * </pre>
60  *
61  * The area in this classes is computing by finding an accurate approximation
62  * to the area integrand using a discrete sine transform fitting \e N equally
63  * spaced points in &sigma;. \e N chosen to ensure full accuracy for
64  * <i>b</i>/\e a &isin; [0.01, 100] or \e f &isin; [&minus;99, 0.99].
65  *
66  * The algorithms are described in
67  * - C. F. F. Karney,
68  * <a href="https://arxiv.org/abs/2208.00492">
69  * Geodesics on an arbitrary ellipsoid of revolution</a>, Aug. 2022;
70  * preprint <a href="https://arxiv.org/abs/2208.00492">
71  * arxiv:2208.00492</a>.
72  * .
73  * See \ref geodellip for the formulation. See the documentation on the
74  * Geodesic class for additional information on the geodesic problems.
75  *
76  * Example of use:
77  * \include example-GeodesicExact.cpp
78  *
79  * <a href="GeodSolve.1.html">GeodSolve</a> is a command-line utility
80  * providing access to the functionality of GeodesicExact and
81  * GeodesicLineExact (via the -E option).
82  **********************************************************************/
83 
85  private:
86  typedef Math::real real;
87  friend class GeodesicLineExact;
88  static const unsigned maxit1_ = 20;
89  unsigned maxit2_;
90  real tiny_, tol0_, tol1_, tol2_, tolb_, xthresh_;
91 
92  enum captype {
93  CAP_NONE = 0U,
94  CAP_E = 1U<<0,
95  // Skip 1U<<1 for compatibility with Geodesic (not required)
96  CAP_D = 1U<<2,
97  CAP_H = 1U<<3,
98  CAP_C4 = 1U<<4,
99  CAP_ALL = 0x1FU,
100  CAP_MASK = CAP_ALL,
101  OUT_ALL = 0x7F80U,
102  OUT_MASK = 0xFF80U, // Includes LONG_UNROLL
103  };
104 
105  static real Astroid(real x, real y);
106 
107  real _a, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2;
108  int _nC4;
109  DST _fft;
110 
111  void Lengths(const EllipticFunction& E,
112  real sig12,
113  real ssig1, real csig1, real dn1,
114  real ssig2, real csig2, real dn2,
115  real cbet1, real cbet2, unsigned outmask,
116  real& s12s, real& m12a, real& m0,
117  real& M12, real& M21) const;
118  real InverseStart(EllipticFunction& E,
119  real sbet1, real cbet1, real dn1,
120  real sbet2, real cbet2, real dn2,
121  real lam12, real slam12, real clam12,
122  real& salp1, real& calp1,
123  real& salp2, real& calp2, real& dnm) const;
124  real Lambda12(real sbet1, real cbet1, real dn1,
125  real sbet2, real cbet2, real dn2,
126  real salp1, real calp1, real slam120, real clam120,
127  real& salp2, real& calp2, real& sig12,
128  real& ssig1, real& csig1, real& ssig2, real& csig2,
129  EllipticFunction& E,
130  real& domg12, bool diffp, real& dlam12) const;
131  real GenInverse(real lat1, real lon1, real lat2, real lon2,
132  unsigned outmask, real& s12,
133  real& salp1, real& calp1, real& salp2, real& calp2,
134  real& m12, real& M12, real& M21, real& S12) const;
135 
136  class I4Integrand {
137  private:
138  real X, tX, tdX, sX, sX1, sXX1, asinhsX, _k2;
139  static real asinhsqrt(real x);
140  static real t(real x);
141  static real td(real x);
142  // static real Dt(real x, real y);
143  real DtX(real y) const;
144  public:
145  I4Integrand(real ep2, real k2);
146  real operator()(real sig) const;
147  };
148 
149  public:
150 
151  /**
152  * Bit masks for what calculations to do. These masks do double duty.
153  * They signify to the GeodesicLineExact::GeodesicLineExact constructor and
154  * to GeodesicExact::Line what capabilities should be included in the
155  * GeodesicLineExact object. They also specify which results to return in
156  * the general routines GeodesicExact::GenDirect and
157  * GeodesicExact::GenInverse routines. GeodesicLineExact::mask is a
158  * duplication of this enum.
159  **********************************************************************/
160  enum mask {
161  /**
162  * No capabilities, no output.
163  * @hideinitializer
164  **********************************************************************/
165  NONE = 0U,
166  /**
167  * Calculate latitude \e lat2. (It's not necessary to include this as a
168  * capability to GeodesicLineExact because this is included by default.)
169  * @hideinitializer
170  **********************************************************************/
171  LATITUDE = 1U<<7 | CAP_NONE,
172  /**
173  * Calculate longitude \e lon2.
174  * @hideinitializer
175  **********************************************************************/
176  LONGITUDE = 1U<<8 | CAP_H,
177  /**
178  * Calculate azimuths \e azi1 and \e azi2. (It's not necessary to
179  * include this as a capability to GeodesicLineExact because this is
180  * included by default.)
181  * @hideinitializer
182  **********************************************************************/
183  AZIMUTH = 1U<<9 | CAP_NONE,
184  /**
185  * Calculate distance \e s12.
186  * @hideinitializer
187  **********************************************************************/
188  DISTANCE = 1U<<10 | CAP_E,
189  /**
190  * A combination of the common capabilities: GeodesicExact::LATITUDE,
191  * GeodesicExact::LONGITUDE, GeodesicExact::AZIMUTH,
192  * GeodesicExact::DISTANCE.
193  * @hideinitializer
194  **********************************************************************/
195  STANDARD = LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
196  /**
197  * Allow distance \e s12 to be used as input in the direct geodesic
198  * problem.
199  * @hideinitializer
200  **********************************************************************/
201  DISTANCE_IN = 1U<<11 | CAP_E,
202  /**
203  * Calculate reduced length \e m12.
204  * @hideinitializer
205  **********************************************************************/
206  REDUCEDLENGTH = 1U<<12 | CAP_D,
207  /**
208  * Calculate geodesic scales \e M12 and \e M21.
209  * @hideinitializer
210  **********************************************************************/
211  GEODESICSCALE = 1U<<13 | CAP_D,
212  /**
213  * Calculate area \e S12.
214  * @hideinitializer
215  **********************************************************************/
216  AREA = 1U<<14 | CAP_C4,
217  /**
218  * Unroll \e lon2 in the direct calculation.
219  * @hideinitializer
220  **********************************************************************/
221  LONG_UNROLL = 1U<<15,
222  /**
223  * All capabilities, calculate everything. (GeodesicExact::LONG_UNROLL
224  * is not included in this mask.)
225  * @hideinitializer
226  **********************************************************************/
227  ALL = OUT_ALL| CAP_ALL,
228  };
229 
230  /** \name Constructor
231  **********************************************************************/
232  ///@{
233  /**
234  * Constructor for an ellipsoid with
235  *
236  * @param[in] a equatorial radius (meters).
237  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
238  * Negative \e f gives a prolate ellipsoid.
239  * @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
240  * positive.
241  **********************************************************************/
242  GeodesicExact(real a, real f);
243  ///@}
244 
245  /** \name Direct geodesic problem specified in terms of distance.
246  **********************************************************************/
247  ///@{
248  /**
249  * Perform the direct geodesic calculation where the length of the geodesic
250  * is specified in terms of distance.
251  *
252  * @param[in] lat1 latitude of point 1 (degrees).
253  * @param[in] lon1 longitude of point 1 (degrees).
254  * @param[in] azi1 azimuth at point 1 (degrees).
255  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
256  * signed.
257  * @param[out] lat2 latitude of point 2 (degrees).
258  * @param[out] lon2 longitude of point 2 (degrees).
259  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
260  * @param[out] m12 reduced length of geodesic (meters).
261  * @param[out] M12 geodesic scale of point 2 relative to point 1
262  * (dimensionless).
263  * @param[out] M21 geodesic scale of point 1 relative to point 2
264  * (dimensionless).
265  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
266  * @return \e a12 arc length of between point 1 and point 2 (degrees).
267  *
268  * \e lat1 should be in the range [&minus;90&deg;, 90&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;]. The values of
382  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
383  * 180&deg;].
384  *
385  * If either point is at a pole, the azimuth is defined by keeping the
386  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
387  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
388  * 180&deg; signifies a geodesic which is not a shortest path. (For a
389  * prolate ellipsoid, an additional condition is necessary for a shortest
390  * path: the longitudinal extent must not exceed of 180&deg;.)
391  *
392  * The following functions are overloaded versions of GeodesicExact::Direct
393  * which omit some of the output parameters.
394  **********************************************************************/
395  void ArcDirect(real lat1, real lon1, real azi1, real a12,
396  real& lat2, real& lon2, real& azi2, real& s12,
397  real& m12, real& M12, real& M21, real& S12)
398  const {
399  GenDirect(lat1, lon1, azi1, true, a12,
400  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
401  REDUCEDLENGTH | GEODESICSCALE | AREA,
402  lat2, lon2, azi2, s12, m12, M12, M21, S12);
403  }
404 
405  /**
406  * See the documentation for GeodesicExact::ArcDirect.
407  **********************************************************************/
408  void ArcDirect(real lat1, real lon1, real azi1, real a12,
409  real& lat2, real& lon2) const {
410  real t;
411  GenDirect(lat1, lon1, azi1, true, a12,
412  LATITUDE | LONGITUDE,
413  lat2, lon2, t, t, t, t, t, t);
414  }
415 
416  /**
417  * See the documentation for GeodesicExact::ArcDirect.
418  **********************************************************************/
419  void ArcDirect(real lat1, real lon1, real azi1, real a12,
420  real& lat2, real& lon2, real& azi2) const {
421  real t;
422  GenDirect(lat1, lon1, azi1, true, a12,
423  LATITUDE | LONGITUDE | AZIMUTH,
424  lat2, lon2, azi2, t, t, t, t, t);
425  }
426 
427  /**
428  * See the documentation for GeodesicExact::ArcDirect.
429  **********************************************************************/
430  void ArcDirect(real lat1, real lon1, real azi1, real a12,
431  real& lat2, real& lon2, real& azi2, real& s12)
432  const {
433  real t;
434  GenDirect(lat1, lon1, azi1, true, a12,
435  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
436  lat2, lon2, azi2, s12, t, t, t, t);
437  }
438 
439  /**
440  * See the documentation for GeodesicExact::ArcDirect.
441  **********************************************************************/
442  void ArcDirect(real lat1, real lon1, real azi1, real a12,
443  real& lat2, real& lon2, real& azi2,
444  real& s12, real& m12) const {
445  real t;
446  GenDirect(lat1, lon1, azi1, true, a12,
447  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
448  REDUCEDLENGTH,
449  lat2, lon2, azi2, s12, m12, t, t, t);
450  }
451 
452  /**
453  * See the documentation for GeodesicExact::ArcDirect.
454  **********************************************************************/
455  void ArcDirect(real lat1, real lon1, real azi1, real a12,
456  real& lat2, real& lon2, real& azi2, real& s12,
457  real& M12, real& M21) const {
458  real t;
459  GenDirect(lat1, lon1, azi1, true, a12,
460  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
461  GEODESICSCALE,
462  lat2, lon2, azi2, s12, t, M12, M21, t);
463  }
464 
465  /**
466  * See the documentation for GeodesicExact::ArcDirect.
467  **********************************************************************/
468  void ArcDirect(real lat1, real lon1, real azi1, real a12,
469  real& lat2, real& lon2, real& azi2, real& s12,
470  real& m12, real& M12, real& M21) const {
471  real t;
472  GenDirect(lat1, lon1, azi1, true, a12,
473  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
474  REDUCEDLENGTH | GEODESICSCALE,
475  lat2, lon2, azi2, s12, m12, M12, M21, t);
476  }
477  ///@}
478 
479  /** \name General version of the direct geodesic solution.
480  **********************************************************************/
481  ///@{
482 
483  /**
484  * The general direct geodesic calculation. GeodesicExact::Direct and
485  * GeodesicExact::ArcDirect are defined in terms of this function.
486  *
487  * @param[in] lat1 latitude of point 1 (degrees).
488  * @param[in] lon1 longitude of point 1 (degrees).
489  * @param[in] azi1 azimuth at point 1 (degrees).
490  * @param[in] arcmode boolean flag determining the meaning of the second
491  * parameter.
492  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
493  * point 1 and point 2 (meters); otherwise it is the arc length between
494  * point 1 and point 2 (degrees); it can be signed.
495  * @param[in] outmask a bitor'ed combination of GeodesicExact::mask values
496  * specifying which of the following parameters should be set.
497  * @param[out] lat2 latitude of point 2 (degrees).
498  * @param[out] lon2 longitude of point 2 (degrees).
499  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
500  * @param[out] s12 distance between point 1 and point 2 (meters).
501  * @param[out] m12 reduced length of geodesic (meters).
502  * @param[out] M12 geodesic scale of point 2 relative to point 1
503  * (dimensionless).
504  * @param[out] M21 geodesic scale of point 1 relative to point 2
505  * (dimensionless).
506  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
507  * @return \e a12 arc length of between point 1 and point 2 (degrees).
508  *
509  * The GeodesicExact::mask values possible for \e outmask are
510  * - \e outmask |= GeodesicExact::LATITUDE for the latitude \e lat2;
511  * - \e outmask |= GeodesicExact::LONGITUDE for the latitude \e lon2;
512  * - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2;
513  * - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12;
514  * - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e
515  * m12;
516  * - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e
517  * M12 and \e M21;
518  * - \e outmask |= GeodesicExact::AREA for the area \e S12;
519  * - \e outmask |= GeodesicExact::ALL for all of the above;
520  * - \e outmask |= GeodesicExact::LONG_UNROLL to unroll \e lon2 instead of
521  * wrapping it into the range [&minus;180&deg;, 180&deg;].
522  * .
523  * The function value \e a12 is always computed and returned and this
524  * equals \e s12_a12 is \e arcmode is true. If \e outmask includes
525  * GeodesicExact::DISTANCE and \e arcmode is false, then \e s12 = \e
526  * s12_a12. It is not necessary to include GeodesicExact::DISTANCE_IN in
527  * \e outmask; this is automatically included is \e arcmode is false.
528  *
529  * With the GeodesicExact::LONG_UNROLL bit set, the quantity \e lon2
530  * &minus; \e lon1 indicates how many times and in what sense the geodesic
531  * encircles the ellipsoid.
532  **********************************************************************/
533  Math::real GenDirect(real lat1, real lon1, real azi1,
534  bool arcmode, real s12_a12, unsigned outmask,
535  real& lat2, real& lon2, real& azi2,
536  real& s12, real& m12, real& M12, real& M21,
537  real& S12) const;
538  ///@}
539 
540  /** \name Inverse geodesic problem.
541  **********************************************************************/
542  ///@{
543  /**
544  * Perform the inverse geodesic calculation.
545  *
546  * @param[in] lat1 latitude of point 1 (degrees).
547  * @param[in] lon1 longitude of point 1 (degrees).
548  * @param[in] lat2 latitude of point 2 (degrees).
549  * @param[in] lon2 longitude of point 2 (degrees).
550  * @param[out] s12 distance between point 1 and point 2 (meters).
551  * @param[out] azi1 azimuth at point 1 (degrees).
552  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
553  * @param[out] m12 reduced length of geodesic (meters).
554  * @param[out] M12 geodesic scale of point 2 relative to point 1
555  * (dimensionless).
556  * @param[out] M21 geodesic scale of point 1 relative to point 2
557  * (dimensionless).
558  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
559  * @return \e a12 arc length of between point 1 and point 2 (degrees).
560  *
561  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
562  * The values of \e azi1 and \e azi2 returned are in the range
563  * [&minus;180&deg;, 180&deg;].
564  *
565  * If either point is at a pole, the azimuth is defined by keeping the
566  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
567  * and taking the limit &epsilon; &rarr; 0+.
568  *
569  * The following functions are overloaded versions of
570  * GeodesicExact::Inverse which omit some of the output parameters. Note,
571  * however, that the arc length is always computed and returned as the
572  * function value.
573  **********************************************************************/
574  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
575  real& s12, real& azi1, real& azi2, real& m12,
576  real& M12, real& M21, real& S12) const {
577  return GenInverse(lat1, lon1, lat2, lon2,
578  DISTANCE | AZIMUTH |
579  REDUCEDLENGTH | GEODESICSCALE | AREA,
580  s12, azi1, azi2, m12, M12, M21, S12);
581  }
582 
583  /**
584  * See the documentation for GeodesicExact::Inverse.
585  **********************************************************************/
586  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
587  real& s12) const {
588  real t;
589  return GenInverse(lat1, lon1, lat2, lon2,
590  DISTANCE,
591  s12, t, t, t, t, t, t);
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& azi1, real& azi2) const {
599  real t;
600  return GenInverse(lat1, lon1, lat2, lon2,
601  AZIMUTH,
602  t, azi1, azi2, 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& s12, real& azi1, real& azi2)
610  const {
611  real t;
612  return GenInverse(lat1, lon1, lat2, lon2,
613  DISTANCE | AZIMUTH,
614  s12, azi1, azi2, t, t, t, t);
615  }
616 
617  /**
618  * See the documentation for GeodesicExact::Inverse.
619  **********************************************************************/
620  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
621  real& s12, real& azi1, real& azi2, real& m12)
622  const {
623  real t;
624  return GenInverse(lat1, lon1, lat2, lon2,
625  DISTANCE | AZIMUTH | REDUCEDLENGTH,
626  s12, azi1, azi2, m12, t, t, t);
627  }
628 
629  /**
630  * See the documentation for GeodesicExact::Inverse.
631  **********************************************************************/
632  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
633  real& s12, real& azi1, real& azi2,
634  real& M12, real& M21) const {
635  real t;
636  return GenInverse(lat1, lon1, lat2, lon2,
637  DISTANCE | AZIMUTH | GEODESICSCALE,
638  s12, azi1, azi2, t, M12, M21, t);
639  }
640 
641  /**
642  * See the documentation for GeodesicExact::Inverse.
643  **********************************************************************/
644  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
645  real& s12, real& azi1, real& azi2, real& m12,
646  real& M12, real& M21) const {
647  real t;
648  return GenInverse(lat1, lon1, lat2, lon2,
649  DISTANCE | AZIMUTH |
650  REDUCEDLENGTH | GEODESICSCALE,
651  s12, azi1, azi2, m12, M12, M21, t);
652  }
653  ///@}
654 
655  /** \name General version of inverse geodesic solution.
656  **********************************************************************/
657  ///@{
658  /**
659  * The general inverse geodesic calculation. GeodesicExact::Inverse is
660  * defined in terms of this function.
661  *
662  * @param[in] lat1 latitude of point 1 (degrees).
663  * @param[in] lon1 longitude of point 1 (degrees).
664  * @param[in] lat2 latitude of point 2 (degrees).
665  * @param[in] lon2 longitude of point 2 (degrees).
666  * @param[in] outmask a bitor'ed combination of GeodesicExact::mask values
667  * specifying which of the following parameters should be set.
668  * @param[out] s12 distance between point 1 and point 2 (meters).
669  * @param[out] azi1 azimuth at point 1 (degrees).
670  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
671  * @param[out] m12 reduced length of geodesic (meters).
672  * @param[out] M12 geodesic scale of point 2 relative to point 1
673  * (dimensionless).
674  * @param[out] M21 geodesic scale of point 1 relative to point 2
675  * (dimensionless).
676  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
677  * @return \e a12 arc length of between point 1 and point 2 (degrees).
678  *
679  * The GeodesicExact::mask values possible for \e outmask are
680  * - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12;
681  * - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2;
682  * - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e
683  * m12;
684  * - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e
685  * M12 and \e M21;
686  * - \e outmask |= GeodesicExact::AREA for the area \e S12;
687  * - \e outmask |= GeodesicExact::ALL for all of the above.
688  * .
689  * The arc length is always computed and returned as the function value.
690  **********************************************************************/
691  Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
692  unsigned outmask,
693  real& s12, real& azi1, real& azi2,
694  real& m12, real& M12, real& M21, real& S12) const;
695  ///@}
696 
697  /** \name Interface to GeodesicLineExact.
698  **********************************************************************/
699  ///@{
700 
701  /**
702  * Set up to compute several points on a single geodesic.
703  *
704  * @param[in] lat1 latitude of point 1 (degrees).
705  * @param[in] lon1 longitude of point 1 (degrees).
706  * @param[in] azi1 azimuth at point 1 (degrees).
707  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
708  * specifying the capabilities the GeodesicLineExact object should
709  * possess, i.e., which quantities can be returned in calls to
710  * GeodesicLineExact::Position.
711  * @return a GeodesicLineExact object.
712  *
713  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
714  *
715  * The GeodesicExact::mask values are
716  * - \e caps |= GeodesicExact::LATITUDE for the latitude \e lat2; this is
717  * added automatically;
718  * - \e caps |= GeodesicExact::LONGITUDE for the latitude \e lon2;
719  * - \e caps |= GeodesicExact::AZIMUTH for the azimuth \e azi2; this is
720  * added automatically;
721  * - \e caps |= GeodesicExact::DISTANCE for the distance \e s12;
722  * - \e caps |= GeodesicExact::REDUCEDLENGTH for the reduced length \e m12;
723  * - \e caps |= GeodesicExact::GEODESICSCALE for the geodesic scales \e M12
724  * and \e M21;
725  * - \e caps |= GeodesicExact::AREA for the area \e S12;
726  * - \e caps |= GeodesicExact::DISTANCE_IN permits the length of the
727  * geodesic to be given in terms of \e s12; without this capability the
728  * length can only be specified in terms of arc length;
729  * - \e caps |= GeodesicExact::ALL for all of the above.
730  * .
731  * The default value of \e caps is GeodesicExact::ALL which turns on all
732  * the capabilities.
733  *
734  * If the point is at a pole, the azimuth is defined by keeping \e lon1
735  * fixed, writing \e lat1 = &plusmn;(90 &minus; &epsilon;), and taking the
736  * limit &epsilon; &rarr; 0+.
737  **********************************************************************/
738  GeodesicLineExact Line(real lat1, real lon1, real azi1,
739  unsigned caps = ALL) const;
740 
741  /**
742  * Define a GeodesicLineExact in terms of the inverse geodesic problem.
743  *
744  * @param[in] lat1 latitude of point 1 (degrees).
745  * @param[in] lon1 longitude of point 1 (degrees).
746  * @param[in] lat2 latitude of point 2 (degrees).
747  * @param[in] lon2 longitude of point 2 (degrees).
748  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
749  * specifying the capabilities the GeodesicLineExact object should
750  * possess, i.e., which quantities can be returned in calls to
751  * GeodesicLineExact::Position.
752  * @return a GeodesicLineExact object.
753  *
754  * This function sets point 3 of the GeodesicLineExact to correspond to
755  * point 2 of the inverse geodesic problem.
756  *
757  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
758  **********************************************************************/
759  GeodesicLineExact InverseLine(real lat1, real lon1, real lat2, real lon2,
760  unsigned caps = ALL) const;
761 
762  /**
763  * Define a GeodesicLineExact in terms of the direct geodesic problem
764  * specified in terms of distance.
765  *
766  * @param[in] lat1 latitude of point 1 (degrees).
767  * @param[in] lon1 longitude of point 1 (degrees).
768  * @param[in] azi1 azimuth at point 1 (degrees).
769  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
770  * negative.
771  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
772  * specifying the capabilities the GeodesicLineExact object should
773  * possess, i.e., which quantities can be returned in calls to
774  * GeodesicLineExact::Position.
775  * @return a GeodesicLineExact object.
776  *
777  * This function sets point 3 of the GeodesicLineExact to correspond to
778  * point 2 of the direct geodesic problem.
779  *
780  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
781  **********************************************************************/
782  GeodesicLineExact DirectLine(real lat1, real lon1, real azi1, real s12,
783  unsigned caps = ALL) const;
784 
785  /**
786  * Define a GeodesicLineExact in terms of the direct geodesic problem
787  * specified in terms of arc length.
788  *
789  * @param[in] lat1 latitude of point 1 (degrees).
790  * @param[in] lon1 longitude of point 1 (degrees).
791  * @param[in] azi1 azimuth at point 1 (degrees).
792  * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
793  * be negative.
794  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
795  * specifying the capabilities the GeodesicLineExact object should
796  * possess, i.e., which quantities can be returned in calls to
797  * GeodesicLineExact::Position.
798  * @return a GeodesicLineExact object.
799  *
800  * This function sets point 3 of the GeodesicLineExact to correspond to
801  * point 2 of the direct geodesic problem.
802  *
803  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
804  **********************************************************************/
805  GeodesicLineExact ArcDirectLine(real lat1, real lon1, real azi1, real a12,
806  unsigned caps = ALL) const;
807 
808  /**
809  * Define a GeodesicLineExact in terms of the direct geodesic problem
810  * specified in terms of either distance or arc length.
811  *
812  * @param[in] lat1 latitude of point 1 (degrees).
813  * @param[in] lon1 longitude of point 1 (degrees).
814  * @param[in] azi1 azimuth at point 1 (degrees).
815  * @param[in] arcmode boolean flag determining the meaning of the \e
816  * s12_a12.
817  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
818  * point 1 and point 2 (meters); otherwise it is the arc length between
819  * point 1 and point 2 (degrees); it can be negative.
820  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
821  * specifying the capabilities the GeodesicLineExact object should
822  * possess, i.e., which quantities can be returned in calls to
823  * GeodesicLineExact::Position.
824  * @return a GeodesicLineExact object.
825  *
826  * This function sets point 3 of the GeodesicLineExact to correspond to
827  * point 2 of the direct geodesic problem.
828  *
829  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
830  **********************************************************************/
831  GeodesicLineExact GenDirectLine(real lat1, real lon1, real azi1,
832  bool arcmode, real s12_a12,
833  unsigned caps = ALL) const;
834  ///@}
835 
836  /** \name Inspector functions.
837  **********************************************************************/
838  ///@{
839 
840  /**
841  * @return \e a the equatorial radius of the ellipsoid (meters). This is
842  * the value used in the constructor.
843  **********************************************************************/
844  Math::real EquatorialRadius() const { return _a; }
845 
846  /**
847  * @return \e f the flattening of the ellipsoid. This is the
848  * value used in the constructor.
849  **********************************************************************/
850  Math::real Flattening() const { return _f; }
851 
852  /**
853  * @return total area of ellipsoid in meters<sup>2</sup>. The area of a
854  * polygon encircling a pole can be found by adding
855  * GeodesicExact::EllipsoidArea()/2 to the sum of \e S12 for each side of
856  * the polygon.
857  **********************************************************************/
859  { return 4 * Math::pi() * _c2; }
860  ///@}
861 
862  /**
863  * A global instantiation of GeodesicExact with the parameters for the
864  * WGS84 ellipsoid.
865  **********************************************************************/
866  static const GeodesicExact& WGS84();
867 
868  };
869 
870 } // namespace GeographicLib
871 
872 #endif // GEOGRAPHICLIB_GEODESICEXACT_HPP
Header for GeographicLib::Constants class.
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:67
Header for GeographicLib::DST class.
Header for GeographicLib::EllipticFunction class.
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
Discrete sine transforms.
Definition: DST.hpp:63
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