GeographicLib  2.1
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  * See \ref geodellip for the formulation. See the documentation on the
67  * Geodesic class for additional information on the geodesic problems.
68  *
69  * Example of use:
70  * \include example-GeodesicExact.cpp
71  *
72  * <a href="GeodSolve.1.html">GeodSolve</a> is a command-line utility
73  * providing access to the functionality of GeodesicExact and
74  * GeodesicLineExact (via the -E option).
75  **********************************************************************/
76 
78  private:
79  typedef Math::real real;
80  friend class GeodesicLineExact;
81  static const unsigned maxit1_ = 20;
82  unsigned maxit2_;
83  real tiny_, tol0_, tol1_, tol2_, tolb_, xthresh_;
84 
85  enum captype {
86  CAP_NONE = 0U,
87  CAP_E = 1U<<0,
88  // Skip 1U<<1 for compatibility with Geodesic (not required)
89  CAP_D = 1U<<2,
90  CAP_H = 1U<<3,
91  CAP_C4 = 1U<<4,
92  CAP_ALL = 0x1FU,
93  CAP_MASK = CAP_ALL,
94  OUT_ALL = 0x7F80U,
95  OUT_MASK = 0xFF80U, // Includes LONG_UNROLL
96  };
97 
98  static real Astroid(real x, real y);
99 
100  real _a, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2;
101  int _nC4;
102  DST _fft;
103 
104  void Lengths(const EllipticFunction& E,
105  real sig12,
106  real ssig1, real csig1, real dn1,
107  real ssig2, real csig2, real dn2,
108  real cbet1, real cbet2, unsigned outmask,
109  real& s12s, real& m12a, real& m0,
110  real& M12, real& M21) const;
111  real InverseStart(EllipticFunction& E,
112  real sbet1, real cbet1, real dn1,
113  real sbet2, real cbet2, real dn2,
114  real lam12, real slam12, real clam12,
115  real& salp1, real& calp1,
116  real& salp2, real& calp2, real& dnm) const;
117  real Lambda12(real sbet1, real cbet1, real dn1,
118  real sbet2, real cbet2, real dn2,
119  real salp1, real calp1, real slam120, real clam120,
120  real& salp2, real& calp2, real& sig12,
121  real& ssig1, real& csig1, real& ssig2, real& csig2,
122  EllipticFunction& E,
123  real& domg12, bool diffp, real& dlam12) const;
124  real GenInverse(real lat1, real lon1, real lat2, real lon2,
125  unsigned outmask, real& s12,
126  real& salp1, real& calp1, real& salp2, real& calp2,
127  real& m12, real& M12, real& M21, real& S12) const;
128 
129  class I4Integrand {
130  private:
131  real X, tX, tdX, sX, sX1, sXX1, asinhsX, _k2;
132  static real asinhsqrt(real x);
133  static real t(real x);
134  static real td(real x);
135  // static real Dt(real x, real y);
136  real DtX(real y) const;
137  public:
138  I4Integrand(real ep2, real k2);
139  real operator()(real sig) const;
140  };
141 
142  public:
143 
144  /**
145  * Bit masks for what calculations to do. These masks do double duty.
146  * They signify to the GeodesicLineExact::GeodesicLineExact constructor and
147  * to GeodesicExact::Line what capabilities should be included in the
148  * GeodesicLineExact object. They also specify which results to return in
149  * the general routines GeodesicExact::GenDirect and
150  * GeodesicExact::GenInverse routines. GeodesicLineExact::mask is a
151  * duplication of this enum.
152  **********************************************************************/
153  enum mask {
154  /**
155  * No capabilities, no output.
156  * @hideinitializer
157  **********************************************************************/
158  NONE = 0U,
159  /**
160  * Calculate latitude \e lat2. (It's not necessary to include this as a
161  * capability to GeodesicLineExact because this is included by default.)
162  * @hideinitializer
163  **********************************************************************/
164  LATITUDE = 1U<<7 | CAP_NONE,
165  /**
166  * Calculate longitude \e lon2.
167  * @hideinitializer
168  **********************************************************************/
169  LONGITUDE = 1U<<8 | CAP_H,
170  /**
171  * Calculate azimuths \e azi1 and \e azi2. (It's not necessary to
172  * include this as a capability to GeodesicLineExact because this is
173  * included by default.)
174  * @hideinitializer
175  **********************************************************************/
176  AZIMUTH = 1U<<9 | CAP_NONE,
177  /**
178  * Calculate distance \e s12.
179  * @hideinitializer
180  **********************************************************************/
181  DISTANCE = 1U<<10 | CAP_E,
182  /**
183  * A combination of the common capabilities: GeodesicExact::LATITUDE,
184  * GeodesicExact::LONGITUDE, GeodesicExact::AZIMUTH,
185  * GeodesicExact::DISTANCE.
186  * @hideinitializer
187  **********************************************************************/
188  STANDARD = LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
189  /**
190  * Allow distance \e s12 to be used as input in the direct geodesic
191  * problem.
192  * @hideinitializer
193  **********************************************************************/
194  DISTANCE_IN = 1U<<11 | CAP_E,
195  /**
196  * Calculate reduced length \e m12.
197  * @hideinitializer
198  **********************************************************************/
199  REDUCEDLENGTH = 1U<<12 | CAP_D,
200  /**
201  * Calculate geodesic scales \e M12 and \e M21.
202  * @hideinitializer
203  **********************************************************************/
204  GEODESICSCALE = 1U<<13 | CAP_D,
205  /**
206  * Calculate area \e S12.
207  * @hideinitializer
208  **********************************************************************/
209  AREA = 1U<<14 | CAP_C4,
210  /**
211  * Unroll \e lon2 in the direct calculation.
212  * @hideinitializer
213  **********************************************************************/
214  LONG_UNROLL = 1U<<15,
215  /**
216  * All capabilities, calculate everything. (GeodesicExact::LONG_UNROLL
217  * is not included in this mask.)
218  * @hideinitializer
219  **********************************************************************/
220  ALL = OUT_ALL| CAP_ALL,
221  };
222 
223  /** \name Constructor
224  **********************************************************************/
225  ///@{
226  /**
227  * Constructor for an ellipsoid with
228  *
229  * @param[in] a equatorial radius (meters).
230  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
231  * Negative \e f gives a prolate ellipsoid.
232  * @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
233  * positive.
234  **********************************************************************/
235  GeodesicExact(real a, real f);
236  ///@}
237 
238  /** \name Direct geodesic problem specified in terms of distance.
239  **********************************************************************/
240  ///@{
241  /**
242  * Perform the direct geodesic calculation where the length of the geodesic
243  * is specified in terms of distance.
244  *
245  * @param[in] lat1 latitude of point 1 (degrees).
246  * @param[in] lon1 longitude of point 1 (degrees).
247  * @param[in] azi1 azimuth at point 1 (degrees).
248  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
249  * signed.
250  * @param[out] lat2 latitude of point 2 (degrees).
251  * @param[out] lon2 longitude of point 2 (degrees).
252  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
253  * @param[out] m12 reduced length of geodesic (meters).
254  * @param[out] M12 geodesic scale of point 2 relative to point 1
255  * (dimensionless).
256  * @param[out] M21 geodesic scale of point 1 relative to point 2
257  * (dimensionless).
258  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
259  * @return \e a12 arc length of between point 1 and point 2 (degrees).
260  *
261  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The values of
262  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
263  * 180&deg;].
264  *
265  * If either point is at a pole, the azimuth is defined by keeping the
266  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
267  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
268  * 180&deg; signifies a geodesic which is not a shortest path. (For a
269  * prolate ellipsoid, an additional condition is necessary for a shortest
270  * path: the longitudinal extent must not exceed of 180&deg;.)
271  *
272  * The following functions are overloaded versions of GeodesicExact::Direct
273  * which omit some of the output parameters. Note, however, that the arc
274  * length is always computed and returned as the function value.
275  **********************************************************************/
276  Math::real Direct(real lat1, real lon1, real azi1, real s12,
277  real& lat2, real& lon2, real& azi2,
278  real& m12, real& M12, real& M21, real& S12)
279  const {
280  real t;
281  return GenDirect(lat1, lon1, azi1, false, s12,
282  LATITUDE | LONGITUDE | AZIMUTH |
283  REDUCEDLENGTH | GEODESICSCALE | AREA,
284  lat2, lon2, azi2, t, m12, M12, M21, S12);
285  }
286 
287  /**
288  * See the documentation for GeodesicExact::Direct.
289  **********************************************************************/
290  Math::real Direct(real lat1, real lon1, real azi1, real s12,
291  real& lat2, real& lon2)
292  const {
293  real t;
294  return GenDirect(lat1, lon1, azi1, false, s12,
295  LATITUDE | LONGITUDE,
296  lat2, lon2, t, t, t, t, t, t);
297  }
298 
299  /**
300  * See the documentation for GeodesicExact::Direct.
301  **********************************************************************/
302  Math::real Direct(real lat1, real lon1, real azi1, real s12,
303  real& lat2, real& lon2, real& azi2)
304  const {
305  real t;
306  return GenDirect(lat1, lon1, azi1, false, s12,
307  LATITUDE | LONGITUDE | AZIMUTH,
308  lat2, lon2, azi2, t, t, t, t, t);
309  }
310 
311  /**
312  * See the documentation for GeodesicExact::Direct.
313  **********************************************************************/
314  Math::real Direct(real lat1, real lon1, real azi1, real s12,
315  real& lat2, real& lon2, real& azi2, real& m12)
316  const {
317  real t;
318  return GenDirect(lat1, lon1, azi1, false, s12,
319  LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH,
320  lat2, lon2, azi2, t, m12, t, t, t);
321  }
322 
323  /**
324  * See the documentation for GeodesicExact::Direct.
325  **********************************************************************/
326  Math::real Direct(real lat1, real lon1, real azi1, real s12,
327  real& lat2, real& lon2, real& azi2,
328  real& M12, real& M21)
329  const {
330  real t;
331  return GenDirect(lat1, lon1, azi1, false, s12,
332  LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE,
333  lat2, lon2, azi2, t, t, M12, M21, t);
334  }
335 
336  /**
337  * See the documentation for GeodesicExact::Direct.
338  **********************************************************************/
339  Math::real Direct(real lat1, real lon1, real azi1, real s12,
340  real& lat2, real& lon2, real& azi2,
341  real& m12, real& M12, real& M21)
342  const {
343  real t;
344  return GenDirect(lat1, lon1, azi1, false, s12,
345  LATITUDE | LONGITUDE | AZIMUTH |
346  REDUCEDLENGTH | GEODESICSCALE,
347  lat2, lon2, azi2, t, m12, M12, M21, t);
348  }
349  ///@}
350 
351  /** \name Direct geodesic problem specified in terms of arc length.
352  **********************************************************************/
353  ///@{
354  /**
355  * Perform the direct geodesic calculation where the length of the geodesic
356  * is specified in terms of arc length.
357  *
358  * @param[in] lat1 latitude of point 1 (degrees).
359  * @param[in] lon1 longitude of point 1 (degrees).
360  * @param[in] azi1 azimuth at point 1 (degrees).
361  * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
362  * be signed.
363  * @param[out] lat2 latitude of point 2 (degrees).
364  * @param[out] lon2 longitude of point 2 (degrees).
365  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
366  * @param[out] s12 distance between point 1 and point 2 (meters).
367  * @param[out] m12 reduced length of geodesic (meters).
368  * @param[out] M12 geodesic scale of point 2 relative to point 1
369  * (dimensionless).
370  * @param[out] M21 geodesic scale of point 1 relative to point 2
371  * (dimensionless).
372  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
373  *
374  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The values of
375  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
376  * 180&deg;].
377  *
378  * If either point is at a pole, the azimuth is defined by keeping the
379  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
380  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
381  * 180&deg; signifies a geodesic which is not a shortest path. (For a
382  * prolate ellipsoid, an additional condition is necessary for a shortest
383  * path: the longitudinal extent must not exceed of 180&deg;.)
384  *
385  * The following functions are overloaded versions of GeodesicExact::Direct
386  * which omit some of the output parameters.
387  **********************************************************************/
388  void ArcDirect(real lat1, real lon1, real azi1, real a12,
389  real& lat2, real& lon2, real& azi2, real& s12,
390  real& m12, real& M12, real& M21, real& S12)
391  const {
392  GenDirect(lat1, lon1, azi1, true, a12,
393  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
394  REDUCEDLENGTH | GEODESICSCALE | AREA,
395  lat2, lon2, azi2, s12, m12, M12, M21, S12);
396  }
397 
398  /**
399  * See the documentation for GeodesicExact::ArcDirect.
400  **********************************************************************/
401  void ArcDirect(real lat1, real lon1, real azi1, real a12,
402  real& lat2, real& lon2) const {
403  real t;
404  GenDirect(lat1, lon1, azi1, true, a12,
405  LATITUDE | LONGITUDE,
406  lat2, lon2, t, t, t, t, t, t);
407  }
408 
409  /**
410  * See the documentation for GeodesicExact::ArcDirect.
411  **********************************************************************/
412  void ArcDirect(real lat1, real lon1, real azi1, real a12,
413  real& lat2, real& lon2, real& azi2) const {
414  real t;
415  GenDirect(lat1, lon1, azi1, true, a12,
416  LATITUDE | LONGITUDE | AZIMUTH,
417  lat2, lon2, azi2, t, t, t, t, t);
418  }
419 
420  /**
421  * See the documentation for GeodesicExact::ArcDirect.
422  **********************************************************************/
423  void ArcDirect(real lat1, real lon1, real azi1, real a12,
424  real& lat2, real& lon2, real& azi2, real& s12)
425  const {
426  real t;
427  GenDirect(lat1, lon1, azi1, true, a12,
428  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
429  lat2, lon2, azi2, s12, t, t, t, t);
430  }
431 
432  /**
433  * See the documentation for GeodesicExact::ArcDirect.
434  **********************************************************************/
435  void ArcDirect(real lat1, real lon1, real azi1, real a12,
436  real& lat2, real& lon2, real& azi2,
437  real& s12, real& m12) const {
438  real t;
439  GenDirect(lat1, lon1, azi1, true, a12,
440  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
441  REDUCEDLENGTH,
442  lat2, lon2, azi2, s12, m12, t, t, t);
443  }
444 
445  /**
446  * See the documentation for GeodesicExact::ArcDirect.
447  **********************************************************************/
448  void ArcDirect(real lat1, real lon1, real azi1, real a12,
449  real& lat2, real& lon2, real& azi2, real& s12,
450  real& M12, real& M21) const {
451  real t;
452  GenDirect(lat1, lon1, azi1, true, a12,
453  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
454  GEODESICSCALE,
455  lat2, lon2, azi2, s12, t, M12, M21, t);
456  }
457 
458  /**
459  * See the documentation for GeodesicExact::ArcDirect.
460  **********************************************************************/
461  void ArcDirect(real lat1, real lon1, real azi1, real a12,
462  real& lat2, real& lon2, real& azi2, real& s12,
463  real& m12, real& M12, real& M21) const {
464  real t;
465  GenDirect(lat1, lon1, azi1, true, a12,
466  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
467  REDUCEDLENGTH | GEODESICSCALE,
468  lat2, lon2, azi2, s12, m12, M12, M21, t);
469  }
470  ///@}
471 
472  /** \name General version of the direct geodesic solution.
473  **********************************************************************/
474  ///@{
475 
476  /**
477  * The general direct geodesic calculation. GeodesicExact::Direct and
478  * GeodesicExact::ArcDirect are defined in terms of this function.
479  *
480  * @param[in] lat1 latitude of point 1 (degrees).
481  * @param[in] lon1 longitude of point 1 (degrees).
482  * @param[in] azi1 azimuth at point 1 (degrees).
483  * @param[in] arcmode boolean flag determining the meaning of the second
484  * parameter.
485  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
486  * point 1 and point 2 (meters); otherwise it is the arc length between
487  * point 1 and point 2 (degrees); it can be signed.
488  * @param[in] outmask a bitor'ed combination of GeodesicExact::mask values
489  * specifying which of the following parameters should be set.
490  * @param[out] lat2 latitude of point 2 (degrees).
491  * @param[out] lon2 longitude of point 2 (degrees).
492  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
493  * @param[out] s12 distance between point 1 and point 2 (meters).
494  * @param[out] m12 reduced length of geodesic (meters).
495  * @param[out] M12 geodesic scale of point 2 relative to point 1
496  * (dimensionless).
497  * @param[out] M21 geodesic scale of point 1 relative to point 2
498  * (dimensionless).
499  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
500  * @return \e a12 arc length of between point 1 and point 2 (degrees).
501  *
502  * The GeodesicExact::mask values possible for \e outmask are
503  * - \e outmask |= GeodesicExact::LATITUDE for the latitude \e lat2;
504  * - \e outmask |= GeodesicExact::LONGITUDE for the latitude \e lon2;
505  * - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2;
506  * - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12;
507  * - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e
508  * m12;
509  * - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e
510  * M12 and \e M21;
511  * - \e outmask |= GeodesicExact::AREA for the area \e S12;
512  * - \e outmask |= GeodesicExact::ALL for all of the above;
513  * - \e outmask |= GeodesicExact::LONG_UNROLL to unroll \e lon2 instead of
514  * wrapping it into the range [&minus;180&deg;, 180&deg;].
515  * .
516  * The function value \e a12 is always computed and returned and this
517  * equals \e s12_a12 is \e arcmode is true. If \e outmask includes
518  * GeodesicExact::DISTANCE and \e arcmode is false, then \e s12 = \e
519  * s12_a12. It is not necessary to include GeodesicExact::DISTANCE_IN in
520  * \e outmask; this is automatically included is \e arcmode is false.
521  *
522  * With the GeodesicExact::LONG_UNROLL bit set, the quantity \e lon2
523  * &minus; \e lon1 indicates how many times and in what sense the geodesic
524  * encircles the ellipsoid.
525  **********************************************************************/
526  Math::real GenDirect(real lat1, real lon1, real azi1,
527  bool arcmode, real s12_a12, unsigned outmask,
528  real& lat2, real& lon2, real& azi2,
529  real& s12, real& m12, real& M12, real& M21,
530  real& S12) const;
531  ///@}
532 
533  /** \name Inverse geodesic problem.
534  **********************************************************************/
535  ///@{
536  /**
537  * Perform the inverse geodesic calculation.
538  *
539  * @param[in] lat1 latitude of point 1 (degrees).
540  * @param[in] lon1 longitude of point 1 (degrees).
541  * @param[in] lat2 latitude of point 2 (degrees).
542  * @param[in] lon2 longitude of point 2 (degrees).
543  * @param[out] s12 distance between point 1 and point 2 (meters).
544  * @param[out] azi1 azimuth at point 1 (degrees).
545  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
546  * @param[out] m12 reduced length of geodesic (meters).
547  * @param[out] M12 geodesic scale of point 2 relative to point 1
548  * (dimensionless).
549  * @param[out] M21 geodesic scale of point 1 relative to point 2
550  * (dimensionless).
551  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
552  * @return \e a12 arc length of between point 1 and point 2 (degrees).
553  *
554  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
555  * The values of \e azi1 and \e azi2 returned are in the range
556  * [&minus;180&deg;, 180&deg;].
557  *
558  * If either point is at a pole, the azimuth is defined by keeping the
559  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
560  * and taking the limit &epsilon; &rarr; 0+.
561  *
562  * The following functions are overloaded versions of
563  * GeodesicExact::Inverse which omit some of the output parameters. Note,
564  * however, that the arc length is always computed and returned as the
565  * function value.
566  **********************************************************************/
567  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
568  real& s12, real& azi1, real& azi2, real& m12,
569  real& M12, real& M21, real& S12) const {
570  return GenInverse(lat1, lon1, lat2, lon2,
571  DISTANCE | AZIMUTH |
572  REDUCEDLENGTH | GEODESICSCALE | AREA,
573  s12, azi1, azi2, m12, M12, M21, S12);
574  }
575 
576  /**
577  * See the documentation for GeodesicExact::Inverse.
578  **********************************************************************/
579  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
580  real& s12) const {
581  real t;
582  return GenInverse(lat1, lon1, lat2, lon2,
583  DISTANCE,
584  s12, t, t, t, t, t, t);
585  }
586 
587  /**
588  * See the documentation for GeodesicExact::Inverse.
589  **********************************************************************/
590  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
591  real& azi1, real& azi2) const {
592  real t;
593  return GenInverse(lat1, lon1, lat2, lon2,
594  AZIMUTH,
595  t, azi1, azi2, t, t, t, t);
596  }
597 
598  /**
599  * See the documentation for GeodesicExact::Inverse.
600  **********************************************************************/
601  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
602  real& s12, real& azi1, real& azi2)
603  const {
604  real t;
605  return GenInverse(lat1, lon1, lat2, lon2,
606  DISTANCE | AZIMUTH,
607  s12, azi1, azi2, t, t, t, t);
608  }
609 
610  /**
611  * See the documentation for GeodesicExact::Inverse.
612  **********************************************************************/
613  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
614  real& s12, real& azi1, real& azi2, real& m12)
615  const {
616  real t;
617  return GenInverse(lat1, lon1, lat2, lon2,
618  DISTANCE | AZIMUTH | REDUCEDLENGTH,
619  s12, azi1, azi2, m12, t, t, t);
620  }
621 
622  /**
623  * See the documentation for GeodesicExact::Inverse.
624  **********************************************************************/
625  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
626  real& s12, real& azi1, real& azi2,
627  real& M12, real& M21) const {
628  real t;
629  return GenInverse(lat1, lon1, lat2, lon2,
630  DISTANCE | AZIMUTH | GEODESICSCALE,
631  s12, azi1, azi2, t, M12, M21, t);
632  }
633 
634  /**
635  * See the documentation for GeodesicExact::Inverse.
636  **********************************************************************/
637  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
638  real& s12, real& azi1, real& azi2, real& m12,
639  real& M12, real& M21) const {
640  real t;
641  return GenInverse(lat1, lon1, lat2, lon2,
642  DISTANCE | AZIMUTH |
643  REDUCEDLENGTH | GEODESICSCALE,
644  s12, azi1, azi2, m12, M12, M21, t);
645  }
646  ///@}
647 
648  /** \name General version of inverse geodesic solution.
649  **********************************************************************/
650  ///@{
651  /**
652  * The general inverse geodesic calculation. GeodesicExact::Inverse is
653  * defined in terms of this function.
654  *
655  * @param[in] lat1 latitude of point 1 (degrees).
656  * @param[in] lon1 longitude of point 1 (degrees).
657  * @param[in] lat2 latitude of point 2 (degrees).
658  * @param[in] lon2 longitude of point 2 (degrees).
659  * @param[in] outmask a bitor'ed combination of GeodesicExact::mask values
660  * specifying which of the following parameters should be set.
661  * @param[out] s12 distance between point 1 and point 2 (meters).
662  * @param[out] azi1 azimuth at point 1 (degrees).
663  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
664  * @param[out] m12 reduced length of geodesic (meters).
665  * @param[out] M12 geodesic scale of point 2 relative to point 1
666  * (dimensionless).
667  * @param[out] M21 geodesic scale of point 1 relative to point 2
668  * (dimensionless).
669  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
670  * @return \e a12 arc length of between point 1 and point 2 (degrees).
671  *
672  * The GeodesicExact::mask values possible for \e outmask are
673  * - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12;
674  * - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2;
675  * - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e
676  * m12;
677  * - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e
678  * M12 and \e M21;
679  * - \e outmask |= GeodesicExact::AREA for the area \e S12;
680  * - \e outmask |= GeodesicExact::ALL for all of the above.
681  * .
682  * The arc length is always computed and returned as the function value.
683  **********************************************************************/
684  Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
685  unsigned outmask,
686  real& s12, real& azi1, real& azi2,
687  real& m12, real& M12, real& M21, real& S12) const;
688  ///@}
689 
690  /** \name Interface to GeodesicLineExact.
691  **********************************************************************/
692  ///@{
693 
694  /**
695  * Set up to compute several points on a single geodesic.
696  *
697  * @param[in] lat1 latitude of point 1 (degrees).
698  * @param[in] lon1 longitude of point 1 (degrees).
699  * @param[in] azi1 azimuth at point 1 (degrees).
700  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
701  * specifying the capabilities the GeodesicLineExact object should
702  * possess, i.e., which quantities can be returned in calls to
703  * GeodesicLineExact::Position.
704  * @return a GeodesicLineExact object.
705  *
706  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
707  *
708  * The GeodesicExact::mask values are
709  * - \e caps |= GeodesicExact::LATITUDE for the latitude \e lat2; this is
710  * added automatically;
711  * - \e caps |= GeodesicExact::LONGITUDE for the latitude \e lon2;
712  * - \e caps |= GeodesicExact::AZIMUTH for the azimuth \e azi2; this is
713  * added automatically;
714  * - \e caps |= GeodesicExact::DISTANCE for the distance \e s12;
715  * - \e caps |= GeodesicExact::REDUCEDLENGTH for the reduced length \e m12;
716  * - \e caps |= GeodesicExact::GEODESICSCALE for the geodesic scales \e M12
717  * and \e M21;
718  * - \e caps |= GeodesicExact::AREA for the area \e S12;
719  * - \e caps |= GeodesicExact::DISTANCE_IN permits the length of the
720  * geodesic to be given in terms of \e s12; without this capability the
721  * length can only be specified in terms of arc length;
722  * - \e caps |= GeodesicExact::ALL for all of the above.
723  * .
724  * The default value of \e caps is GeodesicExact::ALL which turns on all
725  * the capabilities.
726  *
727  * If the point is at a pole, the azimuth is defined by keeping \e lon1
728  * fixed, writing \e lat1 = &plusmn;(90 &minus; &epsilon;), and taking the
729  * limit &epsilon; &rarr; 0+.
730  **********************************************************************/
731  GeodesicLineExact Line(real lat1, real lon1, real azi1,
732  unsigned caps = ALL) const;
733 
734  /**
735  * Define a GeodesicLineExact in terms of the inverse geodesic problem.
736  *
737  * @param[in] lat1 latitude of point 1 (degrees).
738  * @param[in] lon1 longitude of point 1 (degrees).
739  * @param[in] lat2 latitude of point 2 (degrees).
740  * @param[in] lon2 longitude of point 2 (degrees).
741  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
742  * specifying the capabilities the GeodesicLineExact object should
743  * possess, i.e., which quantities can be returned in calls to
744  * GeodesicLineExact::Position.
745  * @return a GeodesicLineExact object.
746  *
747  * This function sets point 3 of the GeodesicLineExact to correspond to
748  * point 2 of the inverse geodesic problem.
749  *
750  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
751  **********************************************************************/
752  GeodesicLineExact InverseLine(real lat1, real lon1, real lat2, real lon2,
753  unsigned caps = ALL) const;
754 
755  /**
756  * Define a GeodesicLineExact in terms of the direct geodesic problem
757  * specified in terms of distance.
758  *
759  * @param[in] lat1 latitude of point 1 (degrees).
760  * @param[in] lon1 longitude of point 1 (degrees).
761  * @param[in] azi1 azimuth at point 1 (degrees).
762  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
763  * negative.
764  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
765  * specifying the capabilities the GeodesicLineExact object should
766  * possess, i.e., which quantities can be returned in calls to
767  * GeodesicLineExact::Position.
768  * @return a GeodesicLineExact object.
769  *
770  * This function sets point 3 of the GeodesicLineExact to correspond to
771  * point 2 of the direct geodesic problem.
772  *
773  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
774  **********************************************************************/
775  GeodesicLineExact DirectLine(real lat1, real lon1, real azi1, real s12,
776  unsigned caps = ALL) const;
777 
778  /**
779  * Define a GeodesicLineExact in terms of the direct geodesic problem
780  * specified in terms of arc length.
781  *
782  * @param[in] lat1 latitude of point 1 (degrees).
783  * @param[in] lon1 longitude of point 1 (degrees).
784  * @param[in] azi1 azimuth at point 1 (degrees).
785  * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
786  * be negative.
787  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
788  * specifying the capabilities the GeodesicLineExact object should
789  * possess, i.e., which quantities can be returned in calls to
790  * GeodesicLineExact::Position.
791  * @return a GeodesicLineExact object.
792  *
793  * This function sets point 3 of the GeodesicLineExact to correspond to
794  * point 2 of the direct geodesic problem.
795  *
796  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
797  **********************************************************************/
798  GeodesicLineExact ArcDirectLine(real lat1, real lon1, real azi1, real a12,
799  unsigned caps = ALL) const;
800 
801  /**
802  * Define a GeodesicLineExact in terms of the direct geodesic problem
803  * specified in terms of either distance or arc length.
804  *
805  * @param[in] lat1 latitude of point 1 (degrees).
806  * @param[in] lon1 longitude of point 1 (degrees).
807  * @param[in] azi1 azimuth at point 1 (degrees).
808  * @param[in] arcmode boolean flag determining the meaning of the \e
809  * s12_a12.
810  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
811  * point 1 and point 2 (meters); otherwise it is the arc length between
812  * point 1 and point 2 (degrees); it can be negative.
813  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
814  * specifying the capabilities the GeodesicLineExact object should
815  * possess, i.e., which quantities can be returned in calls to
816  * GeodesicLineExact::Position.
817  * @return a GeodesicLineExact object.
818  *
819  * This function sets point 3 of the GeodesicLineExact to correspond to
820  * point 2 of the direct geodesic problem.
821  *
822  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
823  **********************************************************************/
824  GeodesicLineExact GenDirectLine(real lat1, real lon1, real azi1,
825  bool arcmode, real s12_a12,
826  unsigned caps = ALL) const;
827  ///@}
828 
829  /** \name Inspector functions.
830  **********************************************************************/
831  ///@{
832 
833  /**
834  * @return \e a the equatorial radius of the ellipsoid (meters). This is
835  * the value used in the constructor.
836  **********************************************************************/
837  Math::real EquatorialRadius() const { return _a; }
838 
839  /**
840  * @return \e f the flattening of the ellipsoid. This is the
841  * value used in the constructor.
842  **********************************************************************/
843  Math::real Flattening() const { return _f; }
844 
845  /**
846  * @return total area of ellipsoid in meters<sup>2</sup>. The area of a
847  * polygon encircling a pole can be found by adding
848  * GeodesicExact::EllipsoidArea()/2 to the sum of \e S12 for each side of
849  * the polygon.
850  **********************************************************************/
852  { return 4 * Math::pi() * _c2; }
853  ///@}
854 
855  /**
856  * A global instantiation of GeodesicExact with the parameters for the
857  * WGS84 ellipsoid.
858  **********************************************************************/
859  static const GeodesicExact& WGS84();
860 
861  };
862 
863 } // namespace GeographicLib
864 
865 #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:61
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