GeographicLib  2.1.2
Geodesic.hpp
Go to the documentation of this file.
1 /**
2  * \file Geodesic.hpp
3  * \brief Header for GeographicLib::Geodesic class
4  *
5  * Copyright (c) Charles Karney (2009-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_GEODESIC_HPP)
11 #define GEOGRAPHICLIB_GEODESIC_HPP 1
12 
14 
15 #if !defined(GEOGRAPHICLIB_GEODESIC_ORDER)
16 /**
17  * The order of the expansions used by Geodesic.
18  * GEOGRAPHICLIB_GEODESIC_ORDER can be set to any integer in [3, 8].
19  **********************************************************************/
20 # define GEOGRAPHICLIB_GEODESIC_ORDER \
21  (GEOGRAPHICLIB_PRECISION == 2 ? 6 : \
22  (GEOGRAPHICLIB_PRECISION == 1 ? 3 : \
23  (GEOGRAPHICLIB_PRECISION == 3 ? 7 : 8)))
24 #endif
25 
26 namespace GeographicLib {
27 
28  class GeodesicLine;
29 
30  /**
31  * \brief %Geodesic calculations
32  *
33  * The shortest path between two points on an ellipsoid at (\e lat1, \e lon1)
34  * and (\e lat2, \e lon2) is called the geodesic. Its length is \e s12 and
35  * the geodesic from point 1 to point 2 has azimuths \e azi1 and \e azi2 at
36  * the two end points. (The azimuth is the heading measured clockwise from
37  * north. \e azi2 is the "forward" azimuth, i.e., the heading that takes you
38  * beyond point 2 not back to point 1.) In the figure below, latitude is
39  * labeled &phi;, longitude &lambda; (with &lambda;<sub>12</sub> =
40  * &lambda;<sub>2</sub> &minus; &lambda;<sub>1</sub>), and azimuth &alpha;.
41  *
42  * <img src="https://upload.wikimedia.org/wikipedia/commons/c/cb/Geodesic_problem_on_an_ellipsoid.svg" width=250 alt="spheroidal triangle">
43  *
44  * Given \e lat1, \e lon1, \e azi1, and \e s12, we can determine \e lat2, \e
45  * lon2, and \e azi2. This is the \e direct geodesic problem and its
46  * solution is given by the function Geodesic::Direct. (If \e s12 is
47  * sufficiently large that the geodesic wraps more than halfway around the
48  * earth, there will be another geodesic between the points with a smaller \e
49  * s12.)
50  *
51  * Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi1, \e
52  * azi2, and \e s12. This is the \e inverse geodesic problem, whose solution
53  * is given by Geodesic::Inverse. Usually, the solution to the inverse
54  * problem is unique. In cases where there are multiple solutions (all with
55  * the same \e s12, of course), all the solutions can be easily generated
56  * once a particular solution is provided.
57  *
58  * The standard way of specifying the direct problem is the specify the
59  * distance \e s12 to the second point. However it is sometimes useful
60  * instead to specify the arc length \e a12 (in degrees) on the auxiliary
61  * sphere. This is a mathematical construct used in solving the geodesic
62  * problems. The solution of the direct problem in this form is provided by
63  * Geodesic::ArcDirect. An arc length in excess of 180&deg; indicates that
64  * the geodesic is not a shortest path. In addition, the arc length between
65  * an equatorial crossing and the next extremum of latitude for a geodesic is
66  * 90&deg;.
67  *
68  * This class can also calculate several other quantities related to
69  * geodesics. These are:
70  * - <i>reduced length</i>. If we fix the first point and increase \e azi1
71  * by \e dazi1 (radians), the second point is displaced \e m12 \e dazi1 in
72  * the direction \e azi2 + 90&deg;. The quantity \e m12 is called
73  * the "reduced length" and is symmetric under interchange of the two
74  * points. On a curved surface the reduced length obeys a symmetry
75  * relation, \e m12 + \e m21 = 0. On a flat surface, we have \e m12 = \e
76  * s12. The ratio <i>s12</i>/\e m12 gives the azimuthal scale for an
77  * azimuthal equidistant projection.
78  * - <i>geodesic scale</i>. Consider a reference geodesic and a second
79  * geodesic parallel to this one at point 1 and separated by a small
80  * distance \e dt. The separation of the two geodesics at point 2 is \e
81  * M12 \e dt where \e M12 is called the "geodesic scale". \e M21 is
82  * defined similarly (with the geodesics being parallel at point 2). On a
83  * flat surface, we have \e M12 = \e M21 = 1. The quantity 1/\e M12 gives
84  * the scale of the Cassini-Soldner projection.
85  * - <i>area</i>. The area between the geodesic from point 1 to point 2 and
86  * the equation is represented by \e S12; it is the area, measured
87  * counter-clockwise, of the geodesic quadrilateral with corners
88  * (<i>lat1</i>,<i>lon1</i>), (0,<i>lon1</i>), (0,<i>lon2</i>), and
89  * (<i>lat2</i>,<i>lon2</i>). It can be used to compute the area of any
90  * geodesic polygon.
91  *
92  * Overloaded versions of Geodesic::Direct, Geodesic::ArcDirect, and
93  * Geodesic::Inverse allow these quantities to be returned. In addition
94  * there are general functions Geodesic::GenDirect, and Geodesic::GenInverse
95  * which allow an arbitrary set of results to be computed. The quantities \e
96  * m12, \e M12, \e M21 which all specify the behavior of nearby geodesics
97  * obey addition rules. If points 1, 2, and 3 all lie on a single geodesic,
98  * then the following rules hold:
99  * - \e s13 = \e s12 + \e s23
100  * - \e a13 = \e a12 + \e a23
101  * - \e S13 = \e S12 + \e S23
102  * - \e m13 = \e m12 \e M23 + \e m23 \e M21
103  * - \e M13 = \e M12 \e M23 &minus; (1 &minus; \e M12 \e M21) \e m23 / \e m12
104  * - \e M31 = \e M32 \e M21 &minus; (1 &minus; \e M23 \e M32) \e m12 / \e m23
105  *
106  * Additional functionality is provided by the GeodesicLine class, which
107  * allows a sequence of points along a geodesic to be computed.
108  *
109  * The shortest distance returned by the solution of the inverse problem is
110  * (obviously) uniquely defined. However, in a few special cases there are
111  * multiple azimuths which yield the same shortest distance. Here is a
112  * catalog of those cases:
113  * - \e lat1 = &minus;\e lat2 (with neither point at a pole). If \e azi1 =
114  * \e azi2, the geodesic is unique. Otherwise there are two geodesics and
115  * the second one is obtained by setting [\e azi1, \e azi2] &rarr; [\e
116  * azi2, \e azi1], [\e M12, \e M21] &rarr; [\e M21, \e M12], \e S12 &rarr;
117  * &minus;\e S12. (This occurs when the longitude difference is near
118  * &plusmn;180&deg; for oblate ellipsoids.)
119  * - \e lon2 = \e lon1 &plusmn; 180&deg; (with neither point at a pole). If
120  * \e azi1 = 0&deg; or &plusmn;180&deg;, the geodesic is unique. Otherwise
121  * there are two geodesics and the second one is obtained by setting [\e
122  * azi1, \e azi2] &rarr; [&minus;\e azi1, &minus;\e azi2], \e S12 &rarr;
123  * &minus;\e S12. (This occurs when \e lat2 is near &minus;\e lat1 for
124  * prolate ellipsoids.)
125  * - Points 1 and 2 at opposite poles. There are infinitely many geodesics
126  * which can be generated by setting [\e azi1, \e azi2] &rarr; [\e azi1, \e
127  * azi2] + [\e d, &minus;\e d], for arbitrary \e d. (For spheres, this
128  * prescription applies when points 1 and 2 are antipodal.)
129  * - \e s12 = 0 (coincident points). There are infinitely many geodesics
130  * which can be generated by setting [\e azi1, \e azi2] &rarr;
131  * [\e azi1, \e azi2] + [\e d, \e d], for arbitrary \e d.
132  *
133  * The calculations are accurate to better than 15 nm (15 nanometers) for the
134  * WGS84 ellipsoid. See Sec. 9 of
135  * <a href="https://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a> for
136  * details. The algorithms used by this class are based on series expansions
137  * using the flattening \e f as a small parameter. These are only accurate
138  * for |<i>f</i>| &lt; 0.02; however reasonably accurate results will be
139  * obtained for |<i>f</i>| &lt; 0.2. Here is a table of the approximate
140  * maximum error (expressed as a distance) for an ellipsoid with the same
141  * equatorial radius as the WGS84 ellipsoid and different values of the
142  * flattening.<pre>
143  * |f| error
144  * 0.01 25 nm
145  * 0.02 30 nm
146  * 0.05 10 um
147  * 0.1 1.5 mm
148  * 0.2 300 mm
149  * </pre>
150  * For very eccentric ellipsoids, use GeodesicExact instead.
151  *
152  * The algorithms are described in
153  * - C. F. F. Karney,
154  * <a href="https://doi.org/10.1007/s00190-012-0578-z">
155  * Algorithms for geodesics</a>,
156  * J. Geodesy <b>87</b>, 43--55 (2013);
157  * DOI: <a href="https://doi.org/10.1007/s00190-012-0578-z">
158  * 10.1007/s00190-012-0578-z</a>;
159  * addenda:
160  * <a href="https://geographiclib.sourceforge.io/geod-addenda.html">
161  * geod-addenda.html</a>.
162  * .
163  * For more information on geodesics see \ref geodesic.
164  *
165  * Example of use:
166  * \include example-Geodesic.cpp
167  *
168  * <a href="GeodSolve.1.html">GeodSolve</a> is a command-line utility
169  * providing access to the functionality of Geodesic and GeodesicLine.
170  **********************************************************************/
171 
173  private:
174  typedef Math::real real;
175  friend class GeodesicLine;
176  static const int nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
177  static const int nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
178  static const int nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER;
179  static const int nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
180  static const int nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
181  static const int nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
182  static const int nA3x_ = nA3_;
183  static const int nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
184  static const int nC3x_ = (nC3_ * (nC3_ - 1)) / 2;
185  static const int nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER;
186  static const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
187  // Size for temporary array
188  // nC = max(max(nC1_, nC1p_, nC2_) + 1, max(nC3_, nC4_))
189  static const int nC_ = GEOGRAPHICLIB_GEODESIC_ORDER + 1;
190  static const unsigned maxit1_ = 20;
191  unsigned maxit2_;
192  real tiny_, tol0_, tol1_, tol2_, tolb_, xthresh_;
193 
194  enum captype {
195  CAP_NONE = 0U,
196  CAP_C1 = 1U<<0,
197  CAP_C1p = 1U<<1,
198  CAP_C2 = 1U<<2,
199  CAP_C3 = 1U<<3,
200  CAP_C4 = 1U<<4,
201  CAP_ALL = 0x1FU,
202  CAP_MASK = CAP_ALL,
203  OUT_ALL = 0x7F80U,
204  OUT_MASK = 0xFF80U, // Includes LONG_UNROLL
205  };
206 
207  static real SinCosSeries(bool sinp,
208  real sinx, real cosx, const real c[], int n);
209  static real Astroid(real x, real y);
210 
211  real _a, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2;
212  real _aA3x[nA3x_], _cC3x[nC3x_], _cC4x[nC4x_];
213 
214  void Lengths(real eps, real sig12,
215  real ssig1, real csig1, real dn1,
216  real ssig2, real csig2, real dn2,
217  real cbet1, real cbet2, unsigned outmask,
218  real& s12s, real& m12a, real& m0,
219  real& M12, real& M21, real Ca[]) const;
220  real InverseStart(real sbet1, real cbet1, real dn1,
221  real sbet2, real cbet2, real dn2,
222  real lam12, real slam12, real clam12,
223  real& salp1, real& calp1,
224  real& salp2, real& calp2, real& dnm,
225  real Ca[]) const;
226  real Lambda12(real sbet1, real cbet1, real dn1,
227  real sbet2, real cbet2, real dn2,
228  real salp1, real calp1, real slam120, real clam120,
229  real& salp2, real& calp2, real& sig12,
230  real& ssig1, real& csig1, real& ssig2, real& csig2,
231  real& eps, real& domg12,
232  bool diffp, real& dlam12, real Ca[]) const;
233  real GenInverse(real lat1, real lon1, real lat2, real lon2,
234  unsigned outmask, real& s12,
235  real& salp1, real& calp1, real& salp2, real& calp2,
236  real& m12, real& M12, real& M21, real& S12) const;
237 
238  // These are Maxima generated functions to provide series approximations to
239  // the integrals for the ellipsoidal geodesic.
240  static real A1m1f(real eps);
241  static void C1f(real eps, real c[]);
242  static void C1pf(real eps, real c[]);
243  static real A2m1f(real eps);
244  static void C2f(real eps, real c[]);
245 
246  void A3coeff();
247  real A3f(real eps) const;
248  void C3coeff();
249  void C3f(real eps, real c[]) const;
250  void C4coeff();
251  void C4f(real k2, real c[]) const;
252 
253  public:
254 
255  /**
256  * Bit masks for what calculations to do. These masks do double duty.
257  * They signify to the GeodesicLine constructor and to
258  * Geodesic::Line what capabilities should be included in the GeodesicLine
259  * object. They also specify which results to return in the general
260  * routines Geodesic::GenDirect and Geodesic::GenInverse routines.
261  * GeodesicLine::mask is a duplication of this enum.
262  **********************************************************************/
263  enum mask {
264  /**
265  * No capabilities, no output.
266  * @hideinitializer
267  **********************************************************************/
268  NONE = 0U,
269  /**
270  * Calculate latitude \e lat2. (It's not necessary to include this as a
271  * capability to GeodesicLine because this is included by default.)
272  * @hideinitializer
273  **********************************************************************/
274  LATITUDE = 1U<<7 | CAP_NONE,
275  /**
276  * Calculate longitude \e lon2.
277  * @hideinitializer
278  **********************************************************************/
279  LONGITUDE = 1U<<8 | CAP_C3,
280  /**
281  * Calculate azimuths \e azi1 and \e azi2. (It's not necessary to
282  * include this as a capability to GeodesicLine because this is included
283  * by default.)
284  * @hideinitializer
285  **********************************************************************/
286  AZIMUTH = 1U<<9 | CAP_NONE,
287  /**
288  * Calculate distance \e s12.
289  * @hideinitializer
290  **********************************************************************/
291  DISTANCE = 1U<<10 | CAP_C1,
292  /**
293  * A combination of the common capabilities: Geodesic::LATITUDE,
294  * Geodesic::LONGITUDE, Geodesic::AZIMUTH, Geodesic::DISTANCE.
295  * @hideinitializer
296  **********************************************************************/
297  STANDARD = LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
298  /**
299  * Allow distance \e s12 to be used as input in the direct geodesic
300  * problem.
301  * @hideinitializer
302  **********************************************************************/
303  DISTANCE_IN = 1U<<11 | CAP_C1 | CAP_C1p,
304  /**
305  * Calculate reduced length \e m12.
306  * @hideinitializer
307  **********************************************************************/
308  REDUCEDLENGTH = 1U<<12 | CAP_C1 | CAP_C2,
309  /**
310  * Calculate geodesic scales \e M12 and \e M21.
311  * @hideinitializer
312  **********************************************************************/
313  GEODESICSCALE = 1U<<13 | CAP_C1 | CAP_C2,
314  /**
315  * Calculate area \e S12.
316  * @hideinitializer
317  **********************************************************************/
318  AREA = 1U<<14 | CAP_C4,
319  /**
320  * Unroll \e lon2 in the direct calculation.
321  * @hideinitializer
322  **********************************************************************/
323  LONG_UNROLL = 1U<<15,
324  /**
325  * All capabilities, calculate everything. (Geodesic::LONG_UNROLL is not
326  * included in this mask.)
327  * @hideinitializer
328  **********************************************************************/
329  ALL = OUT_ALL| CAP_ALL,
330  };
331 
332  /** \name Constructor
333  **********************************************************************/
334  ///@{
335  /**
336  * Constructor for an ellipsoid with
337  *
338  * @param[in] a equatorial radius (meters).
339  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
340  * Negative \e f gives a prolate ellipsoid.
341  * @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
342  * positive.
343  **********************************************************************/
344  Geodesic(real a, real f);
345  ///@}
346 
347  /** \name Direct geodesic problem specified in terms of distance.
348  **********************************************************************/
349  ///@{
350  /**
351  * Solve the direct geodesic problem where the length of the geodesic
352  * is specified in terms of distance.
353  *
354  * @param[in] lat1 latitude of point 1 (degrees).
355  * @param[in] lon1 longitude of point 1 (degrees).
356  * @param[in] azi1 azimuth at point 1 (degrees).
357  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
358  * negative.
359  * @param[out] lat2 latitude of point 2 (degrees).
360  * @param[out] lon2 longitude of point 2 (degrees).
361  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
362  * @param[out] m12 reduced length of geodesic (meters).
363  * @param[out] M12 geodesic scale of point 2 relative to point 1
364  * (dimensionless).
365  * @param[out] M21 geodesic scale of point 1 relative to point 2
366  * (dimensionless).
367  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
368  * @return \e a12 arc length of between point 1 and point 2 (degrees).
369  *
370  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The values of
371  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
372  * 180&deg;].
373  *
374  * If either point is at a pole, the azimuth is defined by keeping the
375  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
376  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
377  * 180&deg; signifies a geodesic which is not a shortest path. (For a
378  * prolate ellipsoid, an additional condition is necessary for a shortest
379  * path: the longitudinal extent must not exceed of 180&deg;.)
380  *
381  * The following functions are overloaded versions of Geodesic::Direct
382  * which omit some of the output parameters. Note, however, that the arc
383  * length is always computed and returned as the function value.
384  **********************************************************************/
385  Math::real Direct(real lat1, real lon1, real azi1, real s12,
386  real& lat2, real& lon2, real& azi2,
387  real& m12, real& M12, real& M21, real& S12)
388  const {
389  real t;
390  return GenDirect(lat1, lon1, azi1, false, s12,
391  LATITUDE | LONGITUDE | AZIMUTH |
392  REDUCEDLENGTH | GEODESICSCALE | AREA,
393  lat2, lon2, azi2, t, m12, M12, M21, S12);
394  }
395 
396  /**
397  * See the documentation for Geodesic::Direct.
398  **********************************************************************/
399  Math::real Direct(real lat1, real lon1, real azi1, real s12,
400  real& lat2, real& lon2)
401  const {
402  real t;
403  return GenDirect(lat1, lon1, azi1, false, s12,
404  LATITUDE | LONGITUDE,
405  lat2, lon2, t, t, t, t, t, t);
406  }
407 
408  /**
409  * See the documentation for Geodesic::Direct.
410  **********************************************************************/
411  Math::real Direct(real lat1, real lon1, real azi1, real s12,
412  real& lat2, real& lon2, real& azi2)
413  const {
414  real t;
415  return GenDirect(lat1, lon1, azi1, false, s12,
416  LATITUDE | LONGITUDE | AZIMUTH,
417  lat2, lon2, azi2, t, t, t, t, t);
418  }
419 
420  /**
421  * See the documentation for Geodesic::Direct.
422  **********************************************************************/
423  Math::real Direct(real lat1, real lon1, real azi1, real s12,
424  real& lat2, real& lon2, real& azi2, real& m12)
425  const {
426  real t;
427  return GenDirect(lat1, lon1, azi1, false, s12,
428  LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH,
429  lat2, lon2, azi2, t, m12, t, t, t);
430  }
431 
432  /**
433  * See the documentation for Geodesic::Direct.
434  **********************************************************************/
435  Math::real Direct(real lat1, real lon1, real azi1, real s12,
436  real& lat2, real& lon2, real& azi2,
437  real& M12, real& M21)
438  const {
439  real t;
440  return GenDirect(lat1, lon1, azi1, false, s12,
441  LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE,
442  lat2, lon2, azi2, t, t, M12, M21, t);
443  }
444 
445  /**
446  * See the documentation for Geodesic::Direct.
447  **********************************************************************/
448  Math::real Direct(real lat1, real lon1, real azi1, real s12,
449  real& lat2, real& lon2, real& azi2,
450  real& m12, real& M12, real& M21)
451  const {
452  real t;
453  return GenDirect(lat1, lon1, azi1, false, s12,
454  LATITUDE | LONGITUDE | AZIMUTH |
455  REDUCEDLENGTH | GEODESICSCALE,
456  lat2, lon2, azi2, t, m12, M12, M21, t);
457  }
458  ///@}
459 
460  /** \name Direct geodesic problem specified in terms of arc length.
461  **********************************************************************/
462  ///@{
463  /**
464  * Solve the direct geodesic problem where the length of the geodesic
465  * is specified in terms of arc length.
466  *
467  * @param[in] lat1 latitude of point 1 (degrees).
468  * @param[in] lon1 longitude of point 1 (degrees).
469  * @param[in] azi1 azimuth at point 1 (degrees).
470  * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
471  * be negative.
472  * @param[out] lat2 latitude of point 2 (degrees).
473  * @param[out] lon2 longitude of point 2 (degrees).
474  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
475  * @param[out] s12 distance between point 1 and point 2 (meters).
476  * @param[out] m12 reduced length of geodesic (meters).
477  * @param[out] M12 geodesic scale of point 2 relative to point 1
478  * (dimensionless).
479  * @param[out] M21 geodesic scale of point 1 relative to point 2
480  * (dimensionless).
481  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
482  *
483  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The values of
484  * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
485  * 180&deg;].
486  *
487  * If either point is at a pole, the azimuth is defined by keeping the
488  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
489  * and taking the limit &epsilon; &rarr; 0+. An arc length greater that
490  * 180&deg; signifies a geodesic which is not a shortest path. (For a
491  * prolate ellipsoid, an additional condition is necessary for a shortest
492  * path: the longitudinal extent must not exceed of 180&deg;.)
493  *
494  * The following functions are overloaded versions of Geodesic::Direct
495  * which omit some of the output parameters.
496  **********************************************************************/
497  void ArcDirect(real lat1, real lon1, real azi1, real a12,
498  real& lat2, real& lon2, real& azi2, real& s12,
499  real& m12, real& M12, real& M21, real& S12)
500  const {
501  GenDirect(lat1, lon1, azi1, true, a12,
502  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
503  REDUCEDLENGTH | GEODESICSCALE | AREA,
504  lat2, lon2, azi2, s12, m12, M12, M21, S12);
505  }
506 
507  /**
508  * See the documentation for Geodesic::ArcDirect.
509  **********************************************************************/
510  void ArcDirect(real lat1, real lon1, real azi1, real a12,
511  real& lat2, real& lon2) const {
512  real t;
513  GenDirect(lat1, lon1, azi1, true, a12,
514  LATITUDE | LONGITUDE,
515  lat2, lon2, t, t, t, t, t, t);
516  }
517 
518  /**
519  * See the documentation for Geodesic::ArcDirect.
520  **********************************************************************/
521  void ArcDirect(real lat1, real lon1, real azi1, real a12,
522  real& lat2, real& lon2, real& azi2) const {
523  real t;
524  GenDirect(lat1, lon1, azi1, true, a12,
525  LATITUDE | LONGITUDE | AZIMUTH,
526  lat2, lon2, azi2, t, t, t, t, t);
527  }
528 
529  /**
530  * See the documentation for Geodesic::ArcDirect.
531  **********************************************************************/
532  void ArcDirect(real lat1, real lon1, real azi1, real a12,
533  real& lat2, real& lon2, real& azi2, real& s12)
534  const {
535  real t;
536  GenDirect(lat1, lon1, azi1, true, a12,
537  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
538  lat2, lon2, azi2, s12, t, t, t, t);
539  }
540 
541  /**
542  * See the documentation for Geodesic::ArcDirect.
543  **********************************************************************/
544  void ArcDirect(real lat1, real lon1, real azi1, real a12,
545  real& lat2, real& lon2, real& azi2,
546  real& s12, real& m12) const {
547  real t;
548  GenDirect(lat1, lon1, azi1, true, a12,
549  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
550  REDUCEDLENGTH,
551  lat2, lon2, azi2, s12, m12, t, t, t);
552  }
553 
554  /**
555  * See the documentation for Geodesic::ArcDirect.
556  **********************************************************************/
557  void ArcDirect(real lat1, real lon1, real azi1, real a12,
558  real& lat2, real& lon2, real& azi2, real& s12,
559  real& M12, real& M21) const {
560  real t;
561  GenDirect(lat1, lon1, azi1, true, a12,
562  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
563  GEODESICSCALE,
564  lat2, lon2, azi2, s12, t, M12, M21, t);
565  }
566 
567  /**
568  * See the documentation for Geodesic::ArcDirect.
569  **********************************************************************/
570  void ArcDirect(real lat1, real lon1, real azi1, real a12,
571  real& lat2, real& lon2, real& azi2, real& s12,
572  real& m12, real& M12, real& M21) const {
573  real t;
574  GenDirect(lat1, lon1, azi1, true, a12,
575  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
576  REDUCEDLENGTH | GEODESICSCALE,
577  lat2, lon2, azi2, s12, m12, M12, M21, t);
578  }
579  ///@}
580 
581  /** \name General version of the direct geodesic solution.
582  **********************************************************************/
583  ///@{
584 
585  /**
586  * The general direct geodesic problem. Geodesic::Direct and
587  * Geodesic::ArcDirect are defined in terms of this function.
588  *
589  * @param[in] lat1 latitude of point 1 (degrees).
590  * @param[in] lon1 longitude of point 1 (degrees).
591  * @param[in] azi1 azimuth at point 1 (degrees).
592  * @param[in] arcmode boolean flag determining the meaning of the \e
593  * s12_a12.
594  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
595  * point 1 and point 2 (meters); otherwise it is the arc length between
596  * point 1 and point 2 (degrees); it can be negative.
597  * @param[in] outmask a bitor'ed combination of Geodesic::mask values
598  * specifying which of the following parameters should be set.
599  * @param[out] lat2 latitude of point 2 (degrees).
600  * @param[out] lon2 longitude of point 2 (degrees).
601  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
602  * @param[out] s12 distance between point 1 and point 2 (meters).
603  * @param[out] m12 reduced length of geodesic (meters).
604  * @param[out] M12 geodesic scale of point 2 relative to point 1
605  * (dimensionless).
606  * @param[out] M21 geodesic scale of point 1 relative to point 2
607  * (dimensionless).
608  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
609  * @return \e a12 arc length of between point 1 and point 2 (degrees).
610  *
611  * The Geodesic::mask values possible for \e outmask are
612  * - \e outmask |= Geodesic::LATITUDE for the latitude \e lat2;
613  * - \e outmask |= Geodesic::LONGITUDE for the latitude \e lon2;
614  * - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2;
615  * - \e outmask |= Geodesic::DISTANCE for the distance \e s12;
616  * - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e
617  * m12;
618  * - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e
619  * M12 and \e M21;
620  * - \e outmask |= Geodesic::AREA for the area \e S12;
621  * - \e outmask |= Geodesic::ALL for all of the above;
622  * - \e outmask |= Geodesic::LONG_UNROLL to unroll \e lon2 instead of
623  * wrapping it into the range [&minus;180&deg;, 180&deg;].
624  * .
625  * The function value \e a12 is always computed and returned and this
626  * equals \e s12_a12 is \e arcmode is true. If \e outmask includes
627  * Geodesic::DISTANCE and \e arcmode is false, then \e s12 = \e s12_a12.
628  * It is not necessary to include Geodesic::DISTANCE_IN in \e outmask; this
629  * is automatically included is \e arcmode is false.
630  *
631  * With the Geodesic::LONG_UNROLL bit set, the quantity \e lon2 &minus; \e
632  * lon1 indicates how many times and in what sense the geodesic encircles
633  * the ellipsoid.
634  **********************************************************************/
635  Math::real GenDirect(real lat1, real lon1, real azi1,
636  bool arcmode, real s12_a12, unsigned outmask,
637  real& lat2, real& lon2, real& azi2,
638  real& s12, real& m12, real& M12, real& M21,
639  real& S12) const;
640  ///@}
641 
642  /** \name Inverse geodesic problem.
643  **********************************************************************/
644  ///@{
645  /**
646  * Solve the inverse geodesic problem.
647  *
648  * @param[in] lat1 latitude of point 1 (degrees).
649  * @param[in] lon1 longitude of point 1 (degrees).
650  * @param[in] lat2 latitude of point 2 (degrees).
651  * @param[in] lon2 longitude of point 2 (degrees).
652  * @param[out] s12 distance between point 1 and point 2 (meters).
653  * @param[out] azi1 azimuth at point 1 (degrees).
654  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
655  * @param[out] m12 reduced length of geodesic (meters).
656  * @param[out] M12 geodesic scale of point 2 relative to point 1
657  * (dimensionless).
658  * @param[out] M21 geodesic scale of point 1 relative to point 2
659  * (dimensionless).
660  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
661  * @return \e a12 arc length of between point 1 and point 2 (degrees).
662  *
663  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
664  * The values of \e azi1 and \e azi2 returned are in the range
665  * [&minus;180&deg;, 180&deg;].
666  *
667  * If either point is at a pole, the azimuth is defined by keeping the
668  * longitude fixed, writing \e lat = &plusmn;(90&deg; &minus; &epsilon;),
669  * and taking the limit &epsilon; &rarr; 0+.
670  *
671  * The solution to the inverse problem is found using Newton's method. If
672  * this fails to converge (this is very unlikely in geodetic applications
673  * but does occur for very eccentric ellipsoids), then the bisection method
674  * is used to refine the solution.
675  *
676  * The following functions are overloaded versions of Geodesic::Inverse
677  * which omit some of the output parameters. Note, however, that the arc
678  * length is always computed and returned as the function value.
679  **********************************************************************/
680  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
681  real& s12, real& azi1, real& azi2, real& m12,
682  real& M12, real& M21, real& S12) const {
683  return GenInverse(lat1, lon1, lat2, lon2,
684  DISTANCE | AZIMUTH |
685  REDUCEDLENGTH | GEODESICSCALE | AREA,
686  s12, azi1, azi2, m12, M12, M21, S12);
687  }
688 
689  /**
690  * See the documentation for Geodesic::Inverse.
691  **********************************************************************/
692  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
693  real& s12) const {
694  real t;
695  return GenInverse(lat1, lon1, lat2, lon2,
696  DISTANCE,
697  s12, t, t, t, t, t, t);
698  }
699 
700  /**
701  * See the documentation for Geodesic::Inverse.
702  **********************************************************************/
703  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
704  real& azi1, real& azi2) const {
705  real t;
706  return GenInverse(lat1, lon1, lat2, lon2,
707  AZIMUTH,
708  t, azi1, azi2, t, t, t, t);
709  }
710 
711  /**
712  * See the documentation for Geodesic::Inverse.
713  **********************************************************************/
714  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
715  real& s12, real& azi1, real& azi2)
716  const {
717  real t;
718  return GenInverse(lat1, lon1, lat2, lon2,
719  DISTANCE | AZIMUTH,
720  s12, azi1, azi2, t, t, t, t);
721  }
722 
723  /**
724  * See the documentation for Geodesic::Inverse.
725  **********************************************************************/
726  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
727  real& s12, real& azi1, real& azi2, real& m12)
728  const {
729  real t;
730  return GenInverse(lat1, lon1, lat2, lon2,
731  DISTANCE | AZIMUTH | REDUCEDLENGTH,
732  s12, azi1, azi2, m12, t, t, t);
733  }
734 
735  /**
736  * See the documentation for Geodesic::Inverse.
737  **********************************************************************/
738  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
739  real& s12, real& azi1, real& azi2,
740  real& M12, real& M21) const {
741  real t;
742  return GenInverse(lat1, lon1, lat2, lon2,
743  DISTANCE | AZIMUTH | GEODESICSCALE,
744  s12, azi1, azi2, t, M12, M21, t);
745  }
746 
747  /**
748  * See the documentation for Geodesic::Inverse.
749  **********************************************************************/
750  Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
751  real& s12, real& azi1, real& azi2, real& m12,
752  real& M12, real& M21) const {
753  real t;
754  return GenInverse(lat1, lon1, lat2, lon2,
755  DISTANCE | AZIMUTH |
756  REDUCEDLENGTH | GEODESICSCALE,
757  s12, azi1, azi2, m12, M12, M21, t);
758  }
759  ///@}
760 
761  /** \name General version of inverse geodesic solution.
762  **********************************************************************/
763  ///@{
764  /**
765  * The general inverse geodesic calculation. Geodesic::Inverse is defined
766  * in terms of this function.
767  *
768  * @param[in] lat1 latitude of point 1 (degrees).
769  * @param[in] lon1 longitude of point 1 (degrees).
770  * @param[in] lat2 latitude of point 2 (degrees).
771  * @param[in] lon2 longitude of point 2 (degrees).
772  * @param[in] outmask a bitor'ed combination of Geodesic::mask values
773  * specifying which of the following parameters should be set.
774  * @param[out] s12 distance between point 1 and point 2 (meters).
775  * @param[out] azi1 azimuth at point 1 (degrees).
776  * @param[out] azi2 (forward) azimuth at point 2 (degrees).
777  * @param[out] m12 reduced length of geodesic (meters).
778  * @param[out] M12 geodesic scale of point 2 relative to point 1
779  * (dimensionless).
780  * @param[out] M21 geodesic scale of point 1 relative to point 2
781  * (dimensionless).
782  * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
783  * @return \e a12 arc length of between point 1 and point 2 (degrees).
784  *
785  * The Geodesic::mask values possible for \e outmask are
786  * - \e outmask |= Geodesic::DISTANCE for the distance \e s12;
787  * - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2;
788  * - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e
789  * m12;
790  * - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e
791  * M12 and \e M21;
792  * - \e outmask |= Geodesic::AREA for the area \e S12;
793  * - \e outmask |= Geodesic::ALL for all of the above.
794  * .
795  * The arc length is always computed and returned as the function value.
796  **********************************************************************/
797  Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
798  unsigned outmask,
799  real& s12, real& azi1, real& azi2,
800  real& m12, real& M12, real& M21, real& S12) const;
801  ///@}
802 
803  /** \name Interface to GeodesicLine.
804  **********************************************************************/
805  ///@{
806 
807  /**
808  * Set up to compute several points on a single geodesic.
809  *
810  * @param[in] lat1 latitude of point 1 (degrees).
811  * @param[in] lon1 longitude of point 1 (degrees).
812  * @param[in] azi1 azimuth at point 1 (degrees).
813  * @param[in] caps bitor'ed combination of Geodesic::mask values
814  * specifying the capabilities the GeodesicLine object should possess,
815  * i.e., which quantities can be returned in calls to
816  * GeodesicLine::Position.
817  * @return a GeodesicLine object.
818  *
819  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
820  *
821  * The Geodesic::mask values are
822  * - \e caps |= Geodesic::LATITUDE for the latitude \e lat2; this is
823  * added automatically;
824  * - \e caps |= Geodesic::LONGITUDE for the latitude \e lon2;
825  * - \e caps |= Geodesic::AZIMUTH for the azimuth \e azi2; this is
826  * added automatically;
827  * - \e caps |= Geodesic::DISTANCE for the distance \e s12;
828  * - \e caps |= Geodesic::REDUCEDLENGTH for the reduced length \e m12;
829  * - \e caps |= Geodesic::GEODESICSCALE for the geodesic scales \e M12
830  * and \e M21;
831  * - \e caps |= Geodesic::AREA for the area \e S12;
832  * - \e caps |= Geodesic::DISTANCE_IN permits the length of the
833  * geodesic to be given in terms of \e s12; without this capability the
834  * length can only be specified in terms of arc length;
835  * - \e caps |= Geodesic::ALL for all of the above.
836  * .
837  * The default value of \e caps is Geodesic::ALL.
838  *
839  * If the point is at a pole, the azimuth is defined by keeping \e lon1
840  * fixed, writing \e lat1 = &plusmn;(90 &minus; &epsilon;), and taking the
841  * limit &epsilon; &rarr; 0+.
842  **********************************************************************/
843  GeodesicLine Line(real lat1, real lon1, real azi1, unsigned caps = ALL)
844  const;
845 
846  /**
847  * Define a GeodesicLine in terms of the inverse geodesic problem.
848  *
849  * @param[in] lat1 latitude of point 1 (degrees).
850  * @param[in] lon1 longitude of point 1 (degrees).
851  * @param[in] lat2 latitude of point 2 (degrees).
852  * @param[in] lon2 longitude of point 2 (degrees).
853  * @param[in] caps bitor'ed combination of Geodesic::mask values
854  * specifying the capabilities the GeodesicLine object should possess,
855  * i.e., which quantities can be returned in calls to
856  * GeodesicLine::Position.
857  * @return a GeodesicLine object.
858  *
859  * This function sets point 3 of the GeodesicLine to correspond to point 2
860  * of the inverse geodesic problem.
861  *
862  * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
863  **********************************************************************/
864  GeodesicLine InverseLine(real lat1, real lon1, real lat2, real lon2,
865  unsigned caps = ALL) const;
866 
867  /**
868  * Define a GeodesicLine in terms of the direct geodesic problem specified
869  * in terms of distance.
870  *
871  * @param[in] lat1 latitude of point 1 (degrees).
872  * @param[in] lon1 longitude of point 1 (degrees).
873  * @param[in] azi1 azimuth at point 1 (degrees).
874  * @param[in] s12 distance between point 1 and point 2 (meters); it can be
875  * negative.
876  * @param[in] caps bitor'ed combination of Geodesic::mask values
877  * specifying the capabilities the GeodesicLine object should possess,
878  * i.e., which quantities can be returned in calls to
879  * GeodesicLine::Position.
880  * @return a GeodesicLine object.
881  *
882  * This function sets point 3 of the GeodesicLine to correspond to point 2
883  * of the direct geodesic problem.
884  *
885  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
886  **********************************************************************/
887  GeodesicLine DirectLine(real lat1, real lon1, real azi1, real s12,
888  unsigned caps = ALL) const;
889 
890  /**
891  * Define a GeodesicLine in terms of the direct geodesic problem specified
892  * in terms of arc length.
893  *
894  * @param[in] lat1 latitude of point 1 (degrees).
895  * @param[in] lon1 longitude of point 1 (degrees).
896  * @param[in] azi1 azimuth at point 1 (degrees).
897  * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
898  * be negative.
899  * @param[in] caps bitor'ed combination of Geodesic::mask values
900  * specifying the capabilities the GeodesicLine object should possess,
901  * i.e., which quantities can be returned in calls to
902  * GeodesicLine::Position.
903  * @return a GeodesicLine object.
904  *
905  * This function sets point 3 of the GeodesicLine to correspond to point 2
906  * of the direct geodesic problem.
907  *
908  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
909  **********************************************************************/
910  GeodesicLine ArcDirectLine(real lat1, real lon1, real azi1, real a12,
911  unsigned caps = ALL) const;
912 
913  /**
914  * Define a GeodesicLine in terms of the direct geodesic problem specified
915  * in terms of either distance or arc length.
916  *
917  * @param[in] lat1 latitude of point 1 (degrees).
918  * @param[in] lon1 longitude of point 1 (degrees).
919  * @param[in] azi1 azimuth at point 1 (degrees).
920  * @param[in] arcmode boolean flag determining the meaning of the \e
921  * s12_a12.
922  * @param[in] s12_a12 if \e arcmode is false, this is the distance between
923  * point 1 and point 2 (meters); otherwise it is the arc length between
924  * point 1 and point 2 (degrees); it can be negative.
925  * @param[in] caps bitor'ed combination of Geodesic::mask values
926  * specifying the capabilities the GeodesicLine object should possess,
927  * i.e., which quantities can be returned in calls to
928  * GeodesicLine::Position.
929  * @return a GeodesicLine object.
930  *
931  * This function sets point 3 of the GeodesicLine to correspond to point 2
932  * of the direct geodesic problem.
933  *
934  * \e lat1 should be in the range [&minus;90&deg;, 90&deg;].
935  **********************************************************************/
936  GeodesicLine GenDirectLine(real lat1, real lon1, real azi1,
937  bool arcmode, real s12_a12,
938  unsigned caps = ALL) const;
939  ///@}
940 
941  /** \name Inspector functions.
942  **********************************************************************/
943  ///@{
944 
945  /**
946  * @return \e a the equatorial radius of the ellipsoid (meters). This is
947  * the value used in the constructor.
948  **********************************************************************/
949  Math::real EquatorialRadius() const { return _a; }
950 
951  /**
952  * @return \e f the flattening of the ellipsoid. This is the
953  * value used in the constructor.
954  **********************************************************************/
955  Math::real Flattening() const { return _f; }
956 
957  /**
958  * @return total area of ellipsoid in meters<sup>2</sup>. The area of a
959  * polygon encircling a pole can be found by adding
960  * Geodesic::EllipsoidArea()/2 to the sum of \e S12 for each side of the
961  * polygon.
962  **********************************************************************/
964  { return 4 * Math::pi() * _c2; }
965  ///@}
966 
967  /**
968  * A global instantiation of Geodesic with the parameters for the WGS84
969  * ellipsoid.
970  **********************************************************************/
971  static const Geodesic& WGS84();
972 
973  };
974 
975 } // namespace GeographicLib
976 
977 #endif // GEOGRAPHICLIB_GEODESIC_HPP
Header for GeographicLib::Constants class.
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:67
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
#define GEOGRAPHICLIB_GEODESIC_ORDER
Definition: Geodesic.hpp:20
Geodesic calculations
Definition: Geodesic.hpp:172
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12) const
Definition: Geodesic.hpp:423
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2) const
Definition: Geodesic.hpp:411
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &M12, real &M21) const
Definition: Geodesic.hpp:557
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12, real &M12, real &M21) const
Definition: Geodesic.hpp:750
Math::real Flattening() const
Definition: Geodesic.hpp:955
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2) const
Definition: Geodesic.hpp:510
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &M12, real &M21) const
Definition: Geodesic.hpp:738
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
Definition: Geodesic.hpp:570
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
Definition: Geodesic.hpp:385
Math::real EquatorialRadius() const
Definition: Geodesic.hpp:949
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12) const
Definition: Geodesic.hpp:726
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12) const
Definition: Geodesic.hpp:532
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &M12, real &M21) const
Definition: Geodesic.hpp:435
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2) const
Definition: Geodesic.hpp:714
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2) const
Definition: Geodesic.hpp:399
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
Definition: Geodesic.hpp:497
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &azi1, real &azi2) const
Definition: Geodesic.hpp:703
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12) const
Definition: Geodesic.hpp:692
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2, real &s12, real &m12) const
Definition: Geodesic.hpp:544
Math::real EllipsoidArea() const
Definition: Geodesic.hpp:963
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12, real &M12, real &M21) const
Definition: Geodesic.hpp:448
void ArcDirect(real lat1, real lon1, real azi1, real a12, real &lat2, real &lon2, real &azi2) const
Definition: Geodesic.hpp:521
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
Definition: Geodesic.hpp:680
static T pi()
Definition: Math.hpp:190
Namespace for GeographicLib.
Definition: Accumulator.cpp:12