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