GeographicLib  1.43
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-2015) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9 
10 #if !defined(GEOGRAPHICLIB_GEODESICEXACT_HPP)
11 #define GEOGRAPHICLIB_GEODESICEXACT_HPP 1
12 
15 
16 #if !defined(GEOGRAPHICLIB_GEODESICEXACT_ORDER)
17 /**
18  * The order of the expansions used by GeodesicExact.
19  **********************************************************************/
20 # define GEOGRAPHICLIB_GEODESICEXACT_ORDER 30
21 #endif
22 
23 namespace GeographicLib {
24 
25  class GeodesicLineExact;
26 
27  /**
28  * \brief Exact geodesic calculations
29  *
30  * The equations for geodesics on an ellipsoid can be expressed in terms of
31  * incomplete elliptic integrals. The Geodesic class expands these integrals
32  * in a series in the flattening \e f and this provides an accurate solution
33  * for \e f &isin; [-0.01, 0.01]. The GeodesicExact class computes the
34  * ellitpic integrals directly and so provides a solution which is valid for
35  * all \e f. However, in practice, its use should be limited to about
36  * <i>b</i>/\e a &isin; [0.01, 100] or \e f &isin; [-99, 0.99].
37  *
38  * For the WGS84 ellipsoid, these classes are 2--3 times \e slower than the
39  * series solution and 2--3 times \e less \e accurate (because it's less easy
40  * to control round-off errors with the elliptic integral formulation); i.e.,
41  * the error is about 40 nm (40 nanometers) instead of 15 nm. However the
42  * error in the series solution scales as <i>f</i><sup>7</sup> while the
43  * error in the elliptic integral solution depends weakly on \e f. If the
44  * quarter meridian distance is 10000 km and the ratio <i>b</i>/\e a = 1
45  * &minus; \e f is varied then the approximate maximum error (expressed as a
46  * distance) is <pre>
47  * 1 - f error (nm)
48  * 1/128 387
49  * 1/64 345
50  * 1/32 269
51  * 1/16 210
52  * 1/8 115
53  * 1/4 69
54  * 1/2 36
55  * 1 15
56  * 2 25
57  * 4 96
58  * 8 318
59  * 16 985
60  * 32 2352
61  * 64 6008
62  * 128 19024
63  * </pre>
64  *
65  * The computation of the area in these classes is via a 30th order series.
66  * This gives accurate results for <i>b</i>/\e a &isin; [1/2, 2]; the
67  * accuracy is about 8 decimal digits for <i>b</i>/\e a &isin; [1/4, 4].
68  *
69  * See \ref geodellip for the formulation. See the documentation on the
70  * Geodesic class for additional information on the geodesic problems.
71  *
72  * Example of use:
73  * \include example-GeodesicExact.cpp
74  *
75  * <a href="GeodSolve.1.html">GeodSolve</a> is a command-line utility
76  * providing access to the functionality of GeodesicExact and
77  * GeodesicLineExact (via the -E option).
78  **********************************************************************/
79 
81  private:
82  typedef Math::real real;
83  friend class GeodesicLineExact;
84  static const int nC4_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER;
85  static const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
86  static const unsigned maxit1_ = 20;
87  unsigned maxit2_;
88  real tiny_, tol0_, tol1_, tol2_, tolb_, xthresh_;
89 
90  enum captype {
91  CAP_NONE = 0U,
92  CAP_E = 1U<<0,
93  // Skip 1U<<1 for compatibility with Geodesic (not required)
94  CAP_D = 1U<<2,
95  CAP_H = 1U<<3,
96  CAP_C4 = 1U<<4,
97  CAP_ALL = 0x1FU,
98  CAP_MASK = CAP_ALL,
99  OUT_ALL = 0x7F80U,
100  OUT_MASK = 0xFF80U, // Includes LONG_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 _C4x[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,
114  real& s12s, real& m12a, real& m0,
115  bool scalep, 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,
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,
125  real& salp2, real& calp2, real& sig12,
126  real& ssig1, real& csig1, real& ssig2, real& csig2,
127  EllipticFunction& E,
128  real& omg12, bool diffp, real& dlam12) const;
129 
130  // These are Maxima generated functions to provide series approximations to
131  // the integrals for the area.
132  void C4coeff();
133  void C4f(real k2, real c[]) const;
134  // Large coefficients are split so that lo contains the low 52 bits and hi
135  // the rest. This choice avoids double rounding with doubles and higher
136  // precision types. float coefficients will suffer double rounding;
137  // however the accuracy is already lousy for floats.
138  static Math::real inline reale(long long hi, long long lo) {
139  using std::ldexp;
140  return ldexp(real(hi), 52) + lo;
141  }
142 
143  public:
144 
145  /**
146  * Bit masks for what calculations to do. These masks do double duty.
147  * They signify to the GeodesicLineExact::GeodesicLineExact constructor and
148  * to GeodesicExact::Line what capabilities should be included in the
149  * GeodesicLineExact object. They also specify which results to return in
150  * the general routines GeodesicExact::GenDirect and
151  * GeodesicExact::GenInverse routines. GeodesicLineExact::mask is a
152  * duplication of this enum.
153  **********************************************************************/
154  enum mask {
155  /**
156  * No capabilities, no output.
157  * @hideinitializer
158  **********************************************************************/
159  NONE = 0U,
160  /**
161  * Calculate latitude \e lat2. (It's not necessary to include this as a
162  * capability to GeodesicLineExact because this is included by default.)
163  * @hideinitializer
164  **********************************************************************/
165  LATITUDE = 1U<<7 | CAP_NONE,
166  /**
167  * Calculate longitude \e lon2.
168  * @hideinitializer
169  **********************************************************************/
170  LONGITUDE = 1U<<8 | CAP_H,
171  /**
172  * Calculate azimuths \e azi1 and \e azi2. (It's not necessary to
173  * include this as a capability to GeodesicLineExact because this is
174  * included by default.)
175  * @hideinitializer
176  **********************************************************************/
177  AZIMUTH = 1U<<9 | CAP_NONE,
178  /**
179  * Calculate distance \e s12.
180  * @hideinitializer
181  **********************************************************************/
182  DISTANCE = 1U<<10 | CAP_E,
183  /**
184  * Allow distance \e s12 to be used as input in the direct geodesic
185  * problem.
186  * @hideinitializer
187  **********************************************************************/
188  DISTANCE_IN = 1U<<11 | CAP_E,
189  /**
190  * Calculate reduced length \e m12.
191  * @hideinitializer
192  **********************************************************************/
193  REDUCEDLENGTH = 1U<<12 | CAP_D,
194  /**
195  * Calculate geodesic scales \e M12 and \e M21.
196  * @hideinitializer
197  **********************************************************************/
198  GEODESICSCALE = 1U<<13 | CAP_D,
199  /**
200  * Calculate area \e S12.
201  * @hideinitializer
202  **********************************************************************/
203  AREA = 1U<<14 | CAP_C4,
204  /**
205  * Unroll \e lon2 in the direct calculation. (This flag used to be
206  * called LONG_NOWRAP.)
207  * @hideinitializer
208  **********************************************************************/
209  LONG_UNROLL = 1U<<15,
210  /// \cond SKIP
211  LONG_NOWRAP = LONG_UNROLL,
212  /// \endcond
213  /**
214  * All capabilities, calculate everything. (LONG_UNROLL is not
215  * included in this mask.)
216  * @hideinitializer
217  **********************************************************************/
218  ALL = OUT_ALL| CAP_ALL,
219  };
220 
221  /** \name Constructor
222  **********************************************************************/
223  ///@{
224  /**
225  * Constructor for a ellipsoid with
226  *
227  * @param[in] a equatorial radius (meters).
228  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
229  * Negative \e f gives a prolate ellipsoid. If \e f &gt; 1, set
230  * flattening to 1/\e f.
231  * @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
232  * positive.
233  **********************************************************************/
234  GeodesicExact(real a, real f);
235  ///@}
236 
237  /** \name Direct geodesic problem specified in terms of distance.
238  **********************************************************************/
239  ///@{
240  /**
241  * Perform the direct geodesic calculation where the length of the geodesic
242  * is specified in terms of distance.
243  *
244  * @param[in] lat1 latitude of point 1 (degrees).
245  * @param[in] lon1 longitude of point 1 (degrees).
246  * @param[in] azi1 azimuth at point 1 (degrees).
247  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
248  * signed.
249  * @param[out] lat2 latitude of point 2 (degrees).
250  * @param[out] lon2 longitude of point 2 (degrees).
251  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
252  * @param[out] m12 reduced length of geodesic (meters).
253  * @param[out] M12 geodesic scale of point 2 relative to point 1
254  * (dimensionless).
255  * @param[out] M21 geodesic scale of point 1 relative to point 2
256  * (dimensionless).
257  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
258  * @return \e a12 arc length of between point 1 and point 2 (degrees).
259  *
260  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
261  * azi1 should be in the range [&minus;540&deg;, 540&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;]; \e lon1 and \e
375  * azi1 should be in the range [&minus;540&deg;, 540&deg;). The values of
376  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
377  * 180&deg;).
378  *
379  * If either point is at a pole, the azimuth is defined by keeping the
380  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
381  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
382  * 180&deg; signifies a geodesic which is not a shortest path. (For a
383  * prolate ellipsoid, an additional condition is necessary for a shortest
384  * path: the longitudinal extent must not exceed of 180&deg;.)
385  *
386  * The following functions are overloaded versions of GeodesicExact::Direct
387  * which omit some of the output parameters.
388  **********************************************************************/
389  void ArcDirect(real lat1, real lon1, real azi1, real a12,
390  real& lat2, real& lon2, real& azi2, real& s12,
391  real& m12, real& M12, real& M21, real& S12)
392  const {
393  GenDirect(lat1, lon1, azi1, true, a12,
394  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
395  REDUCEDLENGTH | GEODESICSCALE | AREA,
396  lat2, lon2, azi2, s12, m12, M12, M21, S12);
397  }
398 
399  /**
400  * See the documentation for GeodesicExact::ArcDirect.
401  **********************************************************************/
402  void ArcDirect(real lat1, real lon1, real azi1, real a12,
403  real& lat2, real& lon2) const {
404  real t;
405  GenDirect(lat1, lon1, azi1, true, a12,
406  LATITUDE | LONGITUDE,
407  lat2, lon2, t, t, t, t, t, t);
408  }
409 
410  /**
411  * See the documentation for GeodesicExact::ArcDirect.
412  **********************************************************************/
413  void ArcDirect(real lat1, real lon1, real azi1, real a12,
414  real& lat2, real& lon2, real& azi2) const {
415  real t;
416  GenDirect(lat1, lon1, azi1, true, a12,
417  LATITUDE | LONGITUDE | AZIMUTH,
418  lat2, lon2, azi2, t, t, t, t, t);
419  }
420 
421  /**
422  * See the documentation for GeodesicExact::ArcDirect.
423  **********************************************************************/
424  void ArcDirect(real lat1, real lon1, real azi1, real a12,
425  real& lat2, real& lon2, real& azi2, real& s12)
426  const {
427  real t;
428  GenDirect(lat1, lon1, azi1, true, a12,
429  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
430  lat2, lon2, azi2, s12, t, t, t, t);
431  }
432 
433  /**
434  * See the documentation for GeodesicExact::ArcDirect.
435  **********************************************************************/
436  void ArcDirect(real lat1, real lon1, real azi1, real a12,
437  real& lat2, real& lon2, real& azi2,
438  real& s12, real& m12) const {
439  real t;
440  GenDirect(lat1, lon1, azi1, true, a12,
441  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
442  REDUCEDLENGTH,
443  lat2, lon2, azi2, s12, m12, t, t, t);
444  }
445 
446  /**
447  * See the documentation for GeodesicExact::ArcDirect.
448  **********************************************************************/
449  void ArcDirect(real lat1, real lon1, real azi1, real a12,
450  real& lat2, real& lon2, real& azi2, real& s12,
451  real& M12, real& M21) const {
452  real t;
453  GenDirect(lat1, lon1, azi1, true, a12,
454  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
455  GEODESICSCALE,
456  lat2, lon2, azi2, s12, t, M12, M21, t);
457  }
458 
459  /**
460  * See the documentation for GeodesicExact::ArcDirect.
461  **********************************************************************/
462  void ArcDirect(real lat1, real lon1, real azi1, real a12,
463  real& lat2, real& lon2, real& azi2, real& s12,
464  real& m12, real& M12, real& M21) const {
465  real t;
466  GenDirect(lat1, lon1, azi1, true, a12,
467  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
468  REDUCEDLENGTH | GEODESICSCALE,
469  lat2, lon2, azi2, s12, m12, M12, M21, t);
470  }
471  ///@}
472 
473  /** \name General version of the direct geodesic solution.
474  **********************************************************************/
475  ///@{
476 
477  /**
478  * The general direct geodesic calculation. GeodesicExact::Direct and
479  * GeodesicExact::ArcDirect are defined in terms of this function.
480  *
481  * @param[in] lat1 latitude of point 1 (degrees).
482  * @param[in] lon1 longitude of point 1 (degrees).
483  * @param[in] azi1 azimuth at point 1 (degrees).
484  * @param[in] arcmode boolean flag determining the meaning of the second
485  * parameter.
486  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
487  * point 1 and point 2 (meters); otherwise it is the arc length between
488  * point 1 and point 2 (degrees); it can be signed.
489  * @param[in] outmask a bitor'ed combination of GeodesicExact::mask values
490  * specifying which of the following parameters should be set.
491  * @param[out] lat2 latitude of point 2 (degrees).
492  * @param[out] lon2 longitude of point 2 (degrees).
493  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
494  * @param[out] s12 distance between point 1 and point 2 (meters).
495  * @param[out] m12 reduced length of geodesic (meters).
496  * @param[out] M12 geodesic scale of point 2 relative to point 1
497  * (dimensionless).
498  * @param[out] M21 geodesic scale of point 1 relative to point 2
499  * (dimensionless).
500  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
501  * @return \e a12 arc length of between point 1 and point 2 (degrees).
502  *
503  * The GeodesicExact::mask values possible for \e outmask are
504  * - \e outmask |= GeodesicExact::LATITUDE for the latitude \e lat2;
505  * - \e outmask |= GeodesicExact::LONGITUDE for the latitude \e lon2;
506  * - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2;
507  * - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12;
508  * - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e
509  * m12;
510  * - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e
511  * M12 and \e M21;
512  * - \e outmask |= GeodesicExact::AREA for the area \e S12;
513  * - \e outmask |= GeodesicExact::ALL for all of the above;
514  * - \e outmask |= GeodesicExact::LONG_UNROLL to unroll \e lon2 instead of
515  * wrapping it into the range [&minus;180&deg;, 180&deg;).
516  * .
517  * The function value \e a12 is always computed and returned and this
518  * equals \e s12_a12 is \e arcmode is true. If \e outmask includes
519  * GeodesicExact::DISTANCE and \e arcmode is false, then \e s12 = \e
520  * s12_a12. It is not necessary to include GeodesicExact::DISTANCE_IN in
521  * \e outmask; this is automatically included is \e arcmode is false.
522  *
523  * With the GeodesicExact::LONG_UNROLL bit set, the quantity \e lon2
524  * &minus; \e lon1 indicates how many times and in what sense the geodesic
525  * encircles the ellipsoid. Because \e lon2 might be outside the normal
526  * allowed range for longitudes, [&minus;540&deg;, 540&deg;), be sure to
527  * normalize it with Math::AngNormalize2 before using it in other
528  * GeographicLib calls.
529  **********************************************************************/
530  Math::real GenDirect(real lat1, real lon1, real azi1,
531  bool arcmode, real s12_a12, unsigned outmask,
532  real& lat2, real& lon2, real& azi2,
533  real& s12, real& m12, real& M12, real& M21,
534  real& S12) const;
535  ///@}
536 
537  /** \name Inverse geodesic problem.
538  **********************************************************************/
539  ///@{
540  /**
541  * Perform the inverse geodesic calculation.
542  *
543  * @param[in] lat1 latitude of point 1 (degrees).
544  * @param[in] lon1 longitude of point 1 (degrees).
545  * @param[in] lat2 latitude of point 2 (degrees).
546  * @param[in] lon2 longitude of point 2 (degrees).
547  * @param[out] s12 distance between point 1 and point 2 (meters).
548  * @param[out] azi1 azimuth at point 1 (degrees).
549  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
550  * @param[out] m12 reduced length of geodesic (meters).
551  * @param[out] M12 geodesic scale of point 2 relative to point 1
552  * (dimensionless).
553  * @param[out] M21 geodesic scale of point 1 relative to point 2
554  * (dimensionless).
555  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
556  * @return \e a12 arc length of between point 1 and point 2 (degrees).
557  *
558  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e
559  * lon1 and \e lon2 should be in the range [&minus;540&deg;, 540&deg;).
560  * The values of \e azi1 and \e azi2 returned are in the range
561  * [&minus;180&deg;, 180&deg;).
562  *
563  * If either point is at a pole, the azimuth is defined by keeping the
564  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
565  * and taking the limit &epsilon; &rarr; 0+.
566  *
567  * The following functions are overloaded versions of GeodesicExact::Inverse
568  * which omit some of the output parameters. Note, however, that the arc
569  * length is always computed and returned as the function value.
570  **********************************************************************/
571  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
572  real& s12, real& azi1, real& azi2, real& m12,
573  real& M12, real& M21, real& S12) const {
574  return GenInverse(lat1, lon1, lat2, lon2,
575  DISTANCE | AZIMUTH |
576  REDUCEDLENGTH | GEODESICSCALE | AREA,
577  s12, azi1, azi2, m12, M12, M21, S12);
578  }
579 
580  /**
581  * See the documentation for GeodesicExact::Inverse.
582  **********************************************************************/
583  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
584  real& s12) const {
585  real t;
586  return GenInverse(lat1, lon1, lat2, lon2,
587  DISTANCE,
588  s12, t, t, t, t, t, t);
589  }
590 
591  /**
592  * See the documentation for GeodesicExact::Inverse.
593  **********************************************************************/
594  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
595  real& azi1, real& azi2) const {
596  real t;
597  return GenInverse(lat1, lon1, lat2, lon2,
598  AZIMUTH,
599  t, azi1, azi2, t, t, t, t);
600  }
601 
602  /**
603  * See the documentation for GeodesicExact::Inverse.
604  **********************************************************************/
605  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
606  real& s12, real& azi1, real& azi2)
607  const {
608  real t;
609  return GenInverse(lat1, lon1, lat2, lon2,
610  DISTANCE | AZIMUTH,
611  s12, azi1, azi2, t, t, t, t);
612  }
613 
614  /**
615  * See the documentation for GeodesicExact::Inverse.
616  **********************************************************************/
617  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
618  real& s12, real& azi1, real& azi2, real& m12)
619  const {
620  real t;
621  return GenInverse(lat1, lon1, lat2, lon2,
622  DISTANCE | AZIMUTH | REDUCEDLENGTH,
623  s12, azi1, azi2, m12, t, t, t);
624  }
625 
626  /**
627  * See the documentation for GeodesicExact::Inverse.
628  **********************************************************************/
629  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
630  real& s12, real& azi1, real& azi2,
631  real& M12, real& M21) const {
632  real t;
633  return GenInverse(lat1, lon1, lat2, lon2,
634  DISTANCE | AZIMUTH | GEODESICSCALE,
635  s12, azi1, azi2, t, M12, M21, t);
636  }
637 
638  /**
639  * See the documentation for GeodesicExact::Inverse.
640  **********************************************************************/
641  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
642  real& s12, real& azi1, real& azi2, real& m12,
643  real& M12, real& M21) const {
644  real t;
645  return GenInverse(lat1, lon1, lat2, lon2,
646  DISTANCE | AZIMUTH |
647  REDUCEDLENGTH | GEODESICSCALE,
648  s12, azi1, azi2, m12, M12, M21, t);
649  }
650  ///@}
651 
652  /** \name General version of inverse geodesic solution.
653  **********************************************************************/
654  ///@{
655  /**
656  * The general inverse geodesic calculation. GeodesicExact::Inverse is
657  * defined in terms of this function.
658  *
659  * @param[in] lat1 latitude of point 1 (degrees).
660  * @param[in] lon1 longitude of point 1 (degrees).
661  * @param[in] lat2 latitude of point 2 (degrees).
662  * @param[in] lon2 longitude of point 2 (degrees).
663  * @param[in] outmask a bitor'ed combination of GeodesicExact::mask values
664  * specifying which of the following parameters should be set.
665  * @param[out] s12 distance between point 1 and point 2 (meters).
666  * @param[out] azi1 azimuth at point 1 (degrees).
667  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
668  * @param[out] m12 reduced length of geodesic (meters).
669  * @param[out] M12 geodesic scale of point 2 relative to point 1
670  * (dimensionless).
671  * @param[out] M21 geodesic scale of point 1 relative to point 2
672  * (dimensionless).
673  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
674  * @return \e a12 arc length of between point 1 and point 2 (degrees).
675  *
676  * The GeodesicExact::mask values possible for \e outmask are
677  * - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12;
678  * - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2;
679  * - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e
680  * m12;
681  * - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e
682  * M12 and \e M21;
683  * - \e outmask |= GeodesicExact::AREA for the area \e S12;
684  * - \e outmask |= GeodesicExact::ALL for all of the above.
685  * .
686  * The arc length is always computed and returned as the function value.
687  **********************************************************************/
688  Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
689  unsigned outmask,
690  real& s12, real& azi1, real& azi2,
691  real& m12, real& M12, real& M21, real& S12)
692  const;
693  ///@}
694 
695  /** \name Interface to GeodesicLineExact.
696  **********************************************************************/
697  ///@{
698 
699  /**
700  * Set up to compute several points on a single geodesic.
701  *
702  * @param[in] lat1 latitude of point 1 (degrees).
703  * @param[in] lon1 longitude of point 1 (degrees).
704  * @param[in] azi1 azimuth at point 1 (degrees).
705  * @param[in] caps bitor'ed combination of GeodesicExact::mask values
706  * specifying the capabilities the GeodesicLineExact object should
707  * possess, i.e., which quantities can be returned in calls to
708  * GeodesicLineExact::Position.
709  * @return a GeodesicLineExact object.
710  *
711  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
712  * azi1 should be in the range [&minus;540&deg;, 540&deg;).
713  *
714  * The GeodesicExact::mask values are
715  * - \e caps |= GeodesicExact::LATITUDE for the latitude \e lat2; this is
716  * added automatically;
717  * - \e caps |= GeodesicExact::LONGITUDE for the latitude \e lon2;
718  * - \e caps |= GeodesicExact::AZIMUTH for the azimuth \e azi2; this is
719  * added automatically;
720  * - \e caps |= GeodesicExact::DISTANCE for the distance \e s12;
721  * - \e caps |= GeodesicExact::REDUCEDLENGTH for the reduced length \e m12;
722  * - \e caps |= GeodesicExact::GEODESICSCALE for the geodesic scales \e M12
723  * and \e M21;
724  * - \e caps |= GeodesicExact::AREA for the area \e S12;
725  * - \e caps |= GeodesicExact::DISTANCE_IN permits the length of the
726  * geodesic to be given in terms of \e s12; without this capability the
727  * length can only be specified in terms of arc length;
728  * - \e caps |= GeodesicExact::ALL for all of the above.
729  * .
730  * The default value of \e caps is GeodesicExact::ALL which turns on all
731  * the capabilities.
732  *
733  * If the point is at a pole, the azimuth is defined by keeping \e lon1
734  * fixed, writing \e lat1 = &plusmn;(90 &minus; &epsilon;), and taking the
735  * limit &epsilon; &rarr; 0+.
736  **********************************************************************/
737  GeodesicLineExact Line(real lat1, real lon1, real azi1, unsigned caps = ALL)
738  const;
739 
740  ///@}
741 
742  /** \name Inspector functions.
743  **********************************************************************/
744  ///@{
745 
746  /**
747  * @return \e a the equatorial radius of the ellipsoid (meters). This is
748  * the value used in the constructor.
749  **********************************************************************/
750  Math::real MajorRadius() const { return _a; }
751 
752  /**
753  * @return \e f the flattening of the ellipsoid. This is the
754  * value used in the constructor.
755  **********************************************************************/
756  Math::real Flattening() const { return _f; }
757 
758  /// \cond SKIP
759  /**
760  * <b>DEPRECATED</b>
761  * @return \e r the inverse flattening of the ellipsoid.
762  **********************************************************************/
763  Math::real InverseFlattening() const { return 1/_f; }
764  /// \endcond
765 
766  /**
767  * @return total area of ellipsoid in meters<sup>2</sup>. The area of a
768  * polygon encircling a pole can be found by adding
769  * GeodesicExact::EllipsoidArea()/2 to the sum of \e S12 for each side of
770  * the polygon.
771  **********************************************************************/
773  { return 4 * Math::pi() * _c2; }
774  ///@}
775 
776  /**
777  * A global instantiation of GeodesicExact with the parameters for the WGS84
778  * ellipsoid.
779  **********************************************************************/
780  static const GeodesicExact& WGS84();
781 
782  };
783 
784 } // namespace GeographicLib
785 
786 #endif // GEOGRAPHICLIB_GEODESICEXACT_HPP
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &m12, real &M12, real &M21, real &S12) const
static T pi()
Definition: Math.hpp:214
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:90
Math::real Flattening() const
Math::real EllipsoidArea() const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12, real &M12, real &M21) const
GeographicLib::Math::real real
Definition: GeodSolve.cpp:32
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2) const
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12, real &M12, real &M21) const
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2) const
Elliptic integrals and functions.
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &M12, real &M21) const
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12, real &M12, real &M21, real &S12) const
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12) const
#define GEOGRAPHICLIB_GEODESICEXACT_ORDER
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2) const
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &M12, real &M21) const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &azi1, real &azi2) const
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &m12) const
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &M12, real &M21) const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12) const
Header for GeographicLib::EllipticFunction class.
Exact geodesic calculations.
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2) const
Header for GeographicLib::Constants class.
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12) const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12) const
Math::real MajorRadius() const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2) const
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &m12, real &M12, real &M21) const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12, real &M12, real &M21, real &S12) const