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