00001 /** 00002 * \file GeodesicExact.hpp 00003 * \brief Header for GeographicLib::GeodesicExact class 00004 * 00005 * Copyright (c) Charles Karney (2012-2014) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_GEODESICEXACT_HPP) 00011 #define GEOGRAPHICLIB_GEODESICEXACT_HPP 1 00012 00013 #include <GeographicLib/Constants.hpp> 00014 #include <GeographicLib/EllipticFunction.hpp> 00015 00016 #if !defined(GEOGRAPHICLIB_GEODESICEXACT_ORDER) 00017 /** 00018 * The order of the expansions used by GeodesicExact. 00019 **********************************************************************/ 00020 # define GEOGRAPHICLIB_GEODESICEXACT_ORDER 30 00021 #endif 00022 00023 namespace GeographicLib { 00024 00025 class GeodesicLineExact; 00026 00027 /** 00028 * \brief Exact geodesic calculations 00029 * 00030 * The equations for geodesics on an ellipsoid can be expressed in terms of 00031 * incomplete elliptic integrals. The Geodesic class expands these integrals 00032 * in a series in the flattening \e f and this provides an accurate solution 00033 * for \e f ∈ [-0.01, 0.01]. The GeodesicExact class computes the 00034 * ellitpic integrals directly and so provides a solution which is valid for 00035 * all \e f. However, in practice, its use should be limited to about 00036 * <i>b</i>/\e a ∈ [0.01, 100] or \e f ∈ [-99, 0.99]. 00037 * 00038 * For the WGS84 ellipsoid, these classes are 2--3 times \e slower than the 00039 * series solution and 2--3 times \e less \e accurate (because it's less easy 00040 * to control round-off errors with the elliptic integral formulation); i.e., 00041 * the error is about 40 nm (40 nanometers) instead of 15 nm. However the 00042 * error in the series solution scales as <i>f</i><sup>7</sup> while the 00043 * error in the elliptic integral solution depends weakly on \e f. If the 00044 * quarter meridian distance is 10000 km and the ratio <i>b</i>/\e a = 1 00045 * − \e f is varied then the approximate maximum error (expressed as a 00046 * distance) is <pre> 00047 * 1 - f error (nm) 00048 * 1/128 387 00049 * 1/64 345 00050 * 1/32 269 00051 * 1/16 210 00052 * 1/8 115 00053 * 1/4 69 00054 * 1/2 36 00055 * 1 15 00056 * 2 25 00057 * 4 96 00058 * 8 318 00059 * 16 985 00060 * 32 2352 00061 * 64 6008 00062 * 128 19024 00063 * </pre> 00064 * 00065 * The computation of the area in these classes is via a 30th order series. 00066 * This gives accurate results for <i>b</i>/\e a ∈ [1/2, 2]; the 00067 * accuracy is about 8 decimal digits for <i>b</i>/\e a ∈ [1/4, 4]. 00068 * 00069 * See \ref geodellip for the formulation. See the documentation on the 00070 * Geodesic class for additional information on the geodesic problems. 00071 * 00072 * Example of use: 00073 * \include example-GeodesicExact.cpp 00074 * 00075 * <a href="GeodSolve.1.html">GeodSolve</a> is a command-line utility 00076 * providing access to the functionality of GeodesicExact and 00077 * GeodesicLineExact (via the -E option). 00078 **********************************************************************/ 00079 00080 class GEOGRAPHICLIB_EXPORT GeodesicExact { 00081 private: 00082 typedef Math::real real; 00083 friend class GeodesicLineExact; 00084 static const int nC4_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER; 00085 static const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2; 00086 static const unsigned maxit1_ = 20; 00087 unsigned maxit2_; 00088 real tiny_, tol0_, tol1_, tol2_, tolb_, xthresh_; 00089 00090 enum captype { 00091 CAP_NONE = 0U, 00092 CAP_E = 1U<<0, 00093 // Skip 1U<<1 for compatibility with Geodesic (not required) 00094 CAP_D = 1U<<2, 00095 CAP_H = 1U<<3, 00096 CAP_C4 = 1U<<4, 00097 CAP_ALL = 0x1FU, 00098 OUT_ALL = 0x7F80U, 00099 }; 00100 00101 static real CosSeries(real sinx, real cosx, const real c[], int n); 00102 static inline real AngRound(real x) { 00103 // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57 00104 // for reals = 0.7 pm on the earth if x is an angle in degrees. (This 00105 // is about 1000 times more resolution than we get with angles around 90 00106 // degrees.) We use this to avoid having to deal with near singular 00107 // cases when x is non-zero but tiny (e.g., 1.0e-200). 00108 using std::abs; 00109 const real z = 1/real(16); 00110 GEOGRAPHICLIB_VOLATILE real y = abs(x); 00111 // The compiler mustn't "simplify" z - (z - y) to y 00112 y = y < z ? z - (z - y) : y; 00113 return x < 0 ? -y : y; 00114 } 00115 static inline void SinCosNorm(real& sinx, real& cosx) { 00116 real r = Math::hypot(sinx, cosx); 00117 sinx /= r; 00118 cosx /= r; 00119 } 00120 static real Astroid(real x, real y); 00121 00122 real _a, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2; 00123 real _C4x[nC4x_]; 00124 00125 void Lengths(const EllipticFunction& E, 00126 real sig12, 00127 real ssig1, real csig1, real dn1, 00128 real ssig2, real csig2, real dn2, 00129 real cbet1, real cbet2, 00130 real& s12s, real& m12a, real& m0, 00131 bool scalep, real& M12, real& M21) const; 00132 real InverseStart(EllipticFunction& E, 00133 real sbet1, real cbet1, real dn1, 00134 real sbet2, real cbet2, real dn2, 00135 real lam12, 00136 real& salp1, real& calp1, 00137 real& salp2, real& calp2, real& dnm) const; 00138 real Lambda12(real sbet1, real cbet1, real dn1, 00139 real sbet2, real cbet2, real dn2, 00140 real salp1, real calp1, 00141 real& salp2, real& calp2, real& sig12, 00142 real& ssig1, real& csig1, real& ssig2, real& csig2, 00143 EllipticFunction& E, 00144 real& omg12, bool diffp, real& dlam12) const; 00145 00146 // These are Maxima generated functions to provide series approximations to 00147 // the integrals for the area. 00148 void C4coeff(); 00149 void C4f(real k2, real c[]) const; 00150 // Large coefficients are split so that lo contains the low 52 bits and hi 00151 // the rest. This choice avoids double rounding with doubles and higher 00152 // precision types. float coefficients will suffer double rounding; 00153 // however the accuracy is already lousy for floats. 00154 static Math::real inline reale(long long hi, long long lo) { 00155 using std::ldexp; 00156 return ldexp(real(hi), 52) + lo; 00157 } 00158 static const Math::real* rawC4coeff(); 00159 00160 public: 00161 00162 /** 00163 * Bit masks for what calculations to do. These masks do double duty. 00164 * They signify to the GeodesicLineExact::GeodesicLineExact constructor and 00165 * to GeodesicExact::Line what capabilities should be included in the 00166 * GeodesicLineExact object. They also specify which results to return in 00167 * the general routines GeodesicExact::GenDirect and 00168 * GeodesicExact::GenInverse routines. GeodesicLineExact::mask is a 00169 * duplication of this enum. 00170 **********************************************************************/ 00171 enum mask { 00172 /** 00173 * No capabilities, no output. 00174 * @hideinitializer 00175 **********************************************************************/ 00176 NONE = 0U, 00177 /** 00178 * Calculate latitude \e lat2. (It's not necessary to include this as a 00179 * capability to GeodesicLineExact because this is included by default.) 00180 * @hideinitializer 00181 **********************************************************************/ 00182 LATITUDE = 1U<<7 | CAP_NONE, 00183 /** 00184 * Calculate longitude \e lon2. 00185 * @hideinitializer 00186 **********************************************************************/ 00187 LONGITUDE = 1U<<8 | CAP_H, 00188 /** 00189 * Calculate azimuths \e azi1 and \e azi2. (It's not necessary to 00190 * include this as a capability to GeodesicLineExact because this is 00191 * included by default.) 00192 * @hideinitializer 00193 **********************************************************************/ 00194 AZIMUTH = 1U<<9 | CAP_NONE, 00195 /** 00196 * Calculate distance \e s12. 00197 * @hideinitializer 00198 **********************************************************************/ 00199 DISTANCE = 1U<<10 | CAP_E, 00200 /** 00201 * Allow distance \e s12 to be used as input in the direct geodesic 00202 * problem. 00203 * @hideinitializer 00204 **********************************************************************/ 00205 DISTANCE_IN = 1U<<11 | CAP_E, 00206 /** 00207 * Calculate reduced length \e m12. 00208 * @hideinitializer 00209 **********************************************************************/ 00210 REDUCEDLENGTH = 1U<<12 | CAP_D, 00211 /** 00212 * Calculate geodesic scales \e M12 and \e M21. 00213 * @hideinitializer 00214 **********************************************************************/ 00215 GEODESICSCALE = 1U<<13 | CAP_D, 00216 /** 00217 * Calculate area \e S12. 00218 * @hideinitializer 00219 **********************************************************************/ 00220 AREA = 1U<<14 | CAP_C4, 00221 /** 00222 * All capabilities, calculate everything. 00223 * @hideinitializer 00224 **********************************************************************/ 00225 ALL = OUT_ALL| CAP_ALL, 00226 }; 00227 00228 /** \name Constructor 00229 **********************************************************************/ 00230 ///@{ 00231 /** 00232 * Constructor for a ellipsoid with 00233 * 00234 * @param[in] a equatorial radius (meters). 00235 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00236 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set 00237 * flattening to 1/\e f. 00238 * @exception GeographicErr if \e a or (1 − \e f) \e a is not 00239 * positive. 00240 **********************************************************************/ 00241 GeodesicExact(real a, real f); 00242 ///@} 00243 00244 /** \name Direct geodesic problem specified in terms of distance. 00245 **********************************************************************/ 00246 ///@{ 00247 /** 00248 * Perform the direct geodesic calculation where the length of the geodesic 00249 * is specified in terms of distance. 00250 * 00251 * @param[in] lat1 latitude of point 1 (degrees). 00252 * @param[in] lon1 longitude of point 1 (degrees). 00253 * @param[in] azi1 azimuth at point 1 (degrees). 00254 * @param[in] s12 distance between point 1 and point 2 (meters); it can be 00255 * signed. 00256 * @param[out] lat2 latitude of point 2 (degrees). 00257 * @param[out] lon2 longitude of point 2 (degrees). 00258 * @param[out] azi2 (forward) azimuth at point 2 (degrees). 00259 * @param[out] m12 reduced length of geodesic (meters). 00260 * @param[out] M12 geodesic scale of point 2 relative to point 1 00261 * (dimensionless). 00262 * @param[out] M21 geodesic scale of point 1 relative to point 2 00263 * (dimensionless). 00264 * @param[out] S12 area under the geodesic (meters<sup>2</sup>). 00265 * @return \e a12 arc length of between point 1 and point 2 (degrees). 00266 * 00267 * \e lat1 should be in the range [−90°, 90°]; \e lon1 and \e 00268 * azi1 should be in the range [−540°, 540°). The values of 00269 * \e lon2 and \e azi2 returned are in the range [−180°, 00270 * 180°). 00271 * 00272 * If either point is at a pole, the azimuth is defined by keeping the 00273 * longitude fixed, writing \e lat = ±(90° − ε), 00274 * and taking the limit ε → 0+. An arc length greater that 00275 * 180° signifies a geodesic which is not a shortest path. (For a 00276 * prolate ellipsoid, an additional condition is necessary for a shortest 00277 * path: the longitudinal extent must not exceed of 180°.) 00278 * 00279 * The following functions are overloaded versions of GeodesicExact::Direct 00280 * which omit some of the output parameters. Note, however, that the arc 00281 * length is always computed and returned as the function value. 00282 **********************************************************************/ 00283 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00284 real& lat2, real& lon2, real& azi2, 00285 real& m12, real& M12, real& M21, real& S12) 00286 const { 00287 real t; 00288 return GenDirect(lat1, lon1, azi1, false, s12, 00289 LATITUDE | LONGITUDE | AZIMUTH | 00290 REDUCEDLENGTH | GEODESICSCALE | AREA, 00291 lat2, lon2, azi2, t, m12, M12, M21, S12); 00292 } 00293 00294 /** 00295 * See the documentation for GeodesicExact::Direct. 00296 **********************************************************************/ 00297 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00298 real& lat2, real& lon2) 00299 const { 00300 real t; 00301 return GenDirect(lat1, lon1, azi1, false, s12, 00302 LATITUDE | LONGITUDE, 00303 lat2, lon2, t, t, t, t, t, t); 00304 } 00305 00306 /** 00307 * See the documentation for GeodesicExact::Direct. 00308 **********************************************************************/ 00309 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00310 real& lat2, real& lon2, real& azi2) 00311 const { 00312 real t; 00313 return GenDirect(lat1, lon1, azi1, false, s12, 00314 LATITUDE | LONGITUDE | AZIMUTH, 00315 lat2, lon2, azi2, t, t, t, t, t); 00316 } 00317 00318 /** 00319 * See the documentation for GeodesicExact::Direct. 00320 **********************************************************************/ 00321 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00322 real& lat2, real& lon2, real& azi2, real& m12) 00323 const { 00324 real t; 00325 return GenDirect(lat1, lon1, azi1, false, s12, 00326 LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH, 00327 lat2, lon2, azi2, t, m12, t, t, t); 00328 } 00329 00330 /** 00331 * See the documentation for GeodesicExact::Direct. 00332 **********************************************************************/ 00333 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00334 real& lat2, real& lon2, real& azi2, 00335 real& M12, real& M21) 00336 const { 00337 real t; 00338 return GenDirect(lat1, lon1, azi1, false, s12, 00339 LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE, 00340 lat2, lon2, azi2, t, t, M12, M21, t); 00341 } 00342 00343 /** 00344 * See the documentation for GeodesicExact::Direct. 00345 **********************************************************************/ 00346 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00347 real& lat2, real& lon2, real& azi2, 00348 real& m12, real& M12, real& M21) 00349 const { 00350 real t; 00351 return GenDirect(lat1, lon1, azi1, false, s12, 00352 LATITUDE | LONGITUDE | AZIMUTH | 00353 REDUCEDLENGTH | GEODESICSCALE, 00354 lat2, lon2, azi2, t, m12, M12, M21, t); 00355 } 00356 ///@} 00357 00358 /** \name Direct geodesic problem specified in terms of arc length. 00359 **********************************************************************/ 00360 ///@{ 00361 /** 00362 * Perform the direct geodesic calculation where the length of the geodesic 00363 * is specified in terms of arc length. 00364 * 00365 * @param[in] lat1 latitude of point 1 (degrees). 00366 * @param[in] lon1 longitude of point 1 (degrees). 00367 * @param[in] azi1 azimuth at point 1 (degrees). 00368 * @param[in] a12 arc length between point 1 and point 2 (degrees); it can 00369 * be signed. 00370 * @param[out] lat2 latitude of point 2 (degrees). 00371 * @param[out] lon2 longitude of point 2 (degrees). 00372 * @param[out] azi2 (forward) azimuth at point 2 (degrees). 00373 * @param[out] s12 distance between point 1 and point 2 (meters). 00374 * @param[out] m12 reduced length of geodesic (meters). 00375 * @param[out] M12 geodesic scale of point 2 relative to point 1 00376 * (dimensionless). 00377 * @param[out] M21 geodesic scale of point 1 relative to point 2 00378 * (dimensionless). 00379 * @param[out] S12 area under the geodesic (meters<sup>2</sup>). 00380 * 00381 * \e lat1 should be in the range [−90°, 90°]; \e lon1 and \e 00382 * azi1 should be in the range [−540°, 540°). The values of 00383 * \e lon2 and \e azi2 returned are in the range [−180°, 00384 * 180°). 00385 * 00386 * If either point is at a pole, the azimuth is defined by keeping the 00387 * longitude fixed, writing \e lat = ±(90° − ε), 00388 * and taking the limit ε → 0+. An arc length greater that 00389 * 180° signifies a geodesic which is not a shortest path. (For a 00390 * prolate ellipsoid, an additional condition is necessary for a shortest 00391 * path: the longitudinal extent must not exceed of 180°.) 00392 * 00393 * The following functions are overloaded versions of GeodesicExact::Direct 00394 * which omit some of the output parameters. 00395 **********************************************************************/ 00396 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00397 real& lat2, real& lon2, real& azi2, real& s12, 00398 real& m12, real& M12, real& M21, real& S12) 00399 const { 00400 GenDirect(lat1, lon1, azi1, true, a12, 00401 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE | 00402 REDUCEDLENGTH | GEODESICSCALE | AREA, 00403 lat2, lon2, azi2, s12, m12, M12, M21, S12); 00404 } 00405 00406 /** 00407 * See the documentation for GeodesicExact::ArcDirect. 00408 **********************************************************************/ 00409 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00410 real& lat2, real& lon2) const { 00411 real t; 00412 GenDirect(lat1, lon1, azi1, true, a12, 00413 LATITUDE | LONGITUDE, 00414 lat2, lon2, t, t, t, t, t, t); 00415 } 00416 00417 /** 00418 * See the documentation for GeodesicExact::ArcDirect. 00419 **********************************************************************/ 00420 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00421 real& lat2, real& lon2, real& azi2) const { 00422 real t; 00423 GenDirect(lat1, lon1, azi1, true, a12, 00424 LATITUDE | LONGITUDE | AZIMUTH, 00425 lat2, lon2, azi2, t, t, t, t, t); 00426 } 00427 00428 /** 00429 * See the documentation for GeodesicExact::ArcDirect. 00430 **********************************************************************/ 00431 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00432 real& lat2, real& lon2, real& azi2, real& s12) 00433 const { 00434 real t; 00435 GenDirect(lat1, lon1, azi1, true, a12, 00436 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE, 00437 lat2, lon2, azi2, s12, t, t, t, t); 00438 } 00439 00440 /** 00441 * See the documentation for GeodesicExact::ArcDirect. 00442 **********************************************************************/ 00443 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00444 real& lat2, real& lon2, real& azi2, 00445 real& s12, real& m12) const { 00446 real t; 00447 GenDirect(lat1, lon1, azi1, true, a12, 00448 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE | 00449 REDUCEDLENGTH, 00450 lat2, lon2, azi2, s12, m12, t, t, t); 00451 } 00452 00453 /** 00454 * See the documentation for GeodesicExact::ArcDirect. 00455 **********************************************************************/ 00456 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00457 real& lat2, real& lon2, real& azi2, real& s12, 00458 real& M12, real& M21) const { 00459 real t; 00460 GenDirect(lat1, lon1, azi1, true, a12, 00461 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE | 00462 GEODESICSCALE, 00463 lat2, lon2, azi2, s12, t, M12, M21, t); 00464 } 00465 00466 /** 00467 * See the documentation for GeodesicExact::ArcDirect. 00468 **********************************************************************/ 00469 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00470 real& lat2, real& lon2, real& azi2, real& s12, 00471 real& m12, real& M12, real& M21) const { 00472 real t; 00473 GenDirect(lat1, lon1, azi1, true, a12, 00474 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE | 00475 REDUCEDLENGTH | GEODESICSCALE, 00476 lat2, lon2, azi2, s12, m12, M12, M21, t); 00477 } 00478 ///@} 00479 00480 /** \name General version of the direct geodesic solution. 00481 **********************************************************************/ 00482 ///@{ 00483 00484 /** 00485 * The general direct geodesic calculation. GeodesicExact::Direct and 00486 * GeodesicExact::ArcDirect are defined in terms of this function. 00487 * 00488 * @param[in] lat1 latitude of point 1 (degrees). 00489 * @param[in] lon1 longitude of point 1 (degrees). 00490 * @param[in] azi1 azimuth at point 1 (degrees). 00491 * @param[in] arcmode boolean flag determining the meaning of the second 00492 * parameter. 00493 * @param[in] s12_a12 if \e arcmode is false, this is the distance between 00494 * point 1 and point 2 (meters); otherwise it is the arc length between 00495 * point 1 and point 2 (degrees); it can be signed. 00496 * @param[in] outmask a bitor'ed combination of GeodesicExact::mask values 00497 * specifying which of the following parameters should be set. 00498 * @param[out] lat2 latitude of point 2 (degrees). 00499 * @param[out] lon2 longitude of point 2 (degrees). 00500 * @param[out] azi2 (forward) azimuth at point 2 (degrees). 00501 * @param[out] s12 distance between point 1 and point 2 (meters). 00502 * @param[out] m12 reduced length of geodesic (meters). 00503 * @param[out] M12 geodesic scale of point 2 relative to point 1 00504 * (dimensionless). 00505 * @param[out] M21 geodesic scale of point 1 relative to point 2 00506 * (dimensionless). 00507 * @param[out] S12 area under the geodesic (meters<sup>2</sup>). 00508 * @return \e a12 arc length of between point 1 and point 2 (degrees). 00509 * 00510 * The GeodesicExact::mask values possible for \e outmask are 00511 * - \e outmask |= GeodesicExact::LATITUDE for the latitude \e lat2; 00512 * - \e outmask |= GeodesicExact::LONGITUDE for the latitude \e lon2; 00513 * - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2; 00514 * - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12; 00515 * - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e 00516 * m12; 00517 * - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e 00518 * M12 and \e M21; 00519 * - \e outmask |= GeodesicExact::AREA for the area \e S12; 00520 * - \e outmask |= GeodesicExact::ALL for all of the above. 00521 * . 00522 * The function value \e a12 is always computed and returned and this 00523 * equals \e s12_a12 is \e arcmode is true. If \e outmask includes 00524 * GeodesicExact::DISTANCE and \e arcmode is false, then \e s12 = \e 00525 * s12_a12. It is not necessary to include GeodesicExact::DISTANCE_IN in 00526 * \e outmask; this is automatically included is \e arcmode is false. 00527 **********************************************************************/ 00528 Math::real GenDirect(real lat1, real lon1, real azi1, 00529 bool arcmode, real s12_a12, unsigned outmask, 00530 real& lat2, real& lon2, real& azi2, 00531 real& s12, real& m12, real& M12, real& M21, 00532 real& S12) const; 00533 ///@} 00534 00535 /** \name Inverse geodesic problem. 00536 **********************************************************************/ 00537 ///@{ 00538 /** 00539 * Perform the inverse geodesic calculation. 00540 * 00541 * @param[in] lat1 latitude of point 1 (degrees). 00542 * @param[in] lon1 longitude of point 1 (degrees). 00543 * @param[in] lat2 latitude of point 2 (degrees). 00544 * @param[in] lon2 longitude of point 2 (degrees). 00545 * @param[out] s12 distance between point 1 and point 2 (meters). 00546 * @param[out] azi1 azimuth at point 1 (degrees). 00547 * @param[out] azi2 (forward) azimuth at point 2 (degrees). 00548 * @param[out] m12 reduced length of geodesic (meters). 00549 * @param[out] M12 geodesic scale of point 2 relative to point 1 00550 * (dimensionless). 00551 * @param[out] M21 geodesic scale of point 1 relative to point 2 00552 * (dimensionless). 00553 * @param[out] S12 area under the geodesic (meters<sup>2</sup>). 00554 * @return \e a12 arc length of between point 1 and point 2 (degrees). 00555 * 00556 * \e lat1 and \e lat2 should be in the range [−90°, 90°]; \e 00557 * lon1 and \e lon2 should be in the range [−540°, 540°). 00558 * The values of \e azi1 and \e azi2 returned are in the range 00559 * [−180°, 180°). 00560 * 00561 * If either point is at a pole, the azimuth is defined by keeping the 00562 * longitude fixed, writing \e lat = ±(90° − ε), 00563 * and taking the limit ε → 0+. 00564 * 00565 * The following functions are overloaded versions of GeodesicExact::Inverse 00566 * which omit some of the output parameters. Note, however, that the arc 00567 * length is always computed and returned as the function value. 00568 **********************************************************************/ 00569 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00570 real& s12, real& azi1, real& azi2, real& m12, 00571 real& M12, real& M21, real& S12) const { 00572 return GenInverse(lat1, lon1, lat2, lon2, 00573 DISTANCE | AZIMUTH | 00574 REDUCEDLENGTH | GEODESICSCALE | AREA, 00575 s12, azi1, azi2, m12, M12, M21, S12); 00576 } 00577 00578 /** 00579 * See the documentation for GeodesicExact::Inverse. 00580 **********************************************************************/ 00581 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00582 real& s12) const { 00583 real t; 00584 return GenInverse(lat1, lon1, lat2, lon2, 00585 DISTANCE, 00586 s12, t, t, t, t, t, t); 00587 } 00588 00589 /** 00590 * See the documentation for GeodesicExact::Inverse. 00591 **********************************************************************/ 00592 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00593 real& azi1, real& azi2) const { 00594 real t; 00595 return GenInverse(lat1, lon1, lat2, lon2, 00596 AZIMUTH, 00597 t, azi1, azi2, t, t, t, t); 00598 } 00599 00600 /** 00601 * See the documentation for GeodesicExact::Inverse. 00602 **********************************************************************/ 00603 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00604 real& s12, real& azi1, real& azi2) 00605 const { 00606 real t; 00607 return GenInverse(lat1, lon1, lat2, lon2, 00608 DISTANCE | AZIMUTH, 00609 s12, azi1, azi2, t, t, t, t); 00610 } 00611 00612 /** 00613 * See the documentation for GeodesicExact::Inverse. 00614 **********************************************************************/ 00615 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00616 real& s12, real& azi1, real& azi2, real& m12) 00617 const { 00618 real t; 00619 return GenInverse(lat1, lon1, lat2, lon2, 00620 DISTANCE | AZIMUTH | REDUCEDLENGTH, 00621 s12, azi1, azi2, m12, t, t, t); 00622 } 00623 00624 /** 00625 * See the documentation for GeodesicExact::Inverse. 00626 **********************************************************************/ 00627 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00628 real& s12, real& azi1, real& azi2, 00629 real& M12, real& M21) const { 00630 real t; 00631 return GenInverse(lat1, lon1, lat2, lon2, 00632 DISTANCE | AZIMUTH | GEODESICSCALE, 00633 s12, azi1, azi2, t, M12, M21, t); 00634 } 00635 00636 /** 00637 * See the documentation for GeodesicExact::Inverse. 00638 **********************************************************************/ 00639 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00640 real& s12, real& azi1, real& azi2, real& m12, 00641 real& M12, real& M21) const { 00642 real t; 00643 return GenInverse(lat1, lon1, lat2, lon2, 00644 DISTANCE | AZIMUTH | 00645 REDUCEDLENGTH | GEODESICSCALE, 00646 s12, azi1, azi2, m12, M12, M21, t); 00647 } 00648 ///@} 00649 00650 /** \name General version of inverse geodesic solution. 00651 **********************************************************************/ 00652 ///@{ 00653 /** 00654 * The general inverse geodesic calculation. GeodesicExact::Inverse is 00655 * defined in terms of this function. 00656 * 00657 * @param[in] lat1 latitude of point 1 (degrees). 00658 * @param[in] lon1 longitude of point 1 (degrees). 00659 * @param[in] lat2 latitude of point 2 (degrees). 00660 * @param[in] lon2 longitude of point 2 (degrees). 00661 * @param[in] outmask a bitor'ed combination of GeodesicExact::mask values 00662 * specifying which of the following parameters should be set. 00663 * @param[out] s12 distance between point 1 and point 2 (meters). 00664 * @param[out] azi1 azimuth at point 1 (degrees). 00665 * @param[out] azi2 (forward) azimuth at point 2 (degrees). 00666 * @param[out] m12 reduced length of geodesic (meters). 00667 * @param[out] M12 geodesic scale of point 2 relative to point 1 00668 * (dimensionless). 00669 * @param[out] M21 geodesic scale of point 1 relative to point 2 00670 * (dimensionless). 00671 * @param[out] S12 area under the geodesic (meters<sup>2</sup>). 00672 * @return \e a12 arc length of between point 1 and point 2 (degrees). 00673 * 00674 * The GeodesicExact::mask values possible for \e outmask are 00675 * - \e outmask |= GeodesicExact::DISTANCE for the distance \e s12; 00676 * - \e outmask |= GeodesicExact::AZIMUTH for the latitude \e azi2; 00677 * - \e outmask |= GeodesicExact::REDUCEDLENGTH for the reduced length \e 00678 * m12; 00679 * - \e outmask |= GeodesicExact::GEODESICSCALE for the geodesic scales \e 00680 * M12 and \e M21; 00681 * - \e outmask |= GeodesicExact::AREA for the area \e S12; 00682 * - \e outmask |= GeodesicExact::ALL for all of the above. 00683 * . 00684 * The arc length is always computed and returned as the function value. 00685 **********************************************************************/ 00686 Math::real GenInverse(real lat1, real lon1, real lat2, real lon2, 00687 unsigned outmask, 00688 real& s12, real& azi1, real& azi2, 00689 real& m12, real& M12, real& M21, real& S12) 00690 const; 00691 ///@} 00692 00693 /** \name Interface to GeodesicLineExact. 00694 **********************************************************************/ 00695 ///@{ 00696 00697 /** 00698 * Set up to compute several points on a single geodesic. 00699 * 00700 * @param[in] lat1 latitude of point 1 (degrees). 00701 * @param[in] lon1 longitude of point 1 (degrees). 00702 * @param[in] azi1 azimuth at point 1 (degrees). 00703 * @param[in] caps bitor'ed combination of GeodesicExact::mask values 00704 * specifying the capabilities the GeodesicLineExact object should 00705 * possess, i.e., which quantities can be returned in calls to 00706 * GeodesicLineExact::Position. 00707 * @return a GeodesicLineExact object. 00708 * 00709 * \e lat1 should be in the range [−90°, 90°]; \e lon1 and \e 00710 * azi1 should be in the range [−540°, 540°). 00711 * 00712 * The GeodesicExact::mask values are 00713 * - \e caps |= GeodesicExact::LATITUDE for the latitude \e lat2; this is 00714 * added automatically; 00715 * - \e caps |= GeodesicExact::LONGITUDE for the latitude \e lon2; 00716 * - \e caps |= GeodesicExact::AZIMUTH for the azimuth \e azi2; this is 00717 * added automatically; 00718 * - \e caps |= GeodesicExact::DISTANCE for the distance \e s12; 00719 * - \e caps |= GeodesicExact::REDUCEDLENGTH for the reduced length \e m12; 00720 * - \e caps |= GeodesicExact::GEODESICSCALE for the geodesic scales \e M12 00721 * and \e M21; 00722 * - \e caps |= GeodesicExact::AREA for the area \e S12; 00723 * - \e caps |= GeodesicExact::DISTANCE_IN permits the length of the 00724 * geodesic to be given in terms of \e s12; without this capability the 00725 * length can only be specified in terms of arc length; 00726 * - \e caps |= GeodesicExact::ALL for all of the above. 00727 * . 00728 * The default value of \e caps is GeodesicExact::ALL which turns on all 00729 * the capabilities. 00730 * 00731 * If the point is at a pole, the azimuth is defined by keeping \e lon1 00732 * fixed, writing \e lat1 = ±(90 − ε), and taking the 00733 * limit ε → 0+. 00734 **********************************************************************/ 00735 GeodesicLineExact Line(real lat1, real lon1, real azi1, unsigned caps = ALL) 00736 const; 00737 00738 ///@} 00739 00740 /** \name Inspector functions. 00741 **********************************************************************/ 00742 ///@{ 00743 00744 /** 00745 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00746 * the value used in the constructor. 00747 **********************************************************************/ 00748 Math::real MajorRadius() const { return _a; } 00749 00750 /** 00751 * @return \e f the flattening of the ellipsoid. This is the 00752 * value used in the constructor. 00753 **********************************************************************/ 00754 Math::real Flattening() const { return _f; } 00755 00756 /// \cond SKIP 00757 /** 00758 * <b>DEPRECATED</b> 00759 * @return \e r the inverse flattening of the ellipsoid. 00760 **********************************************************************/ 00761 Math::real InverseFlattening() const { return 1/_f; } 00762 /// \endcond 00763 00764 /** 00765 * @return total area of ellipsoid in meters<sup>2</sup>. The area of a 00766 * polygon encircling a pole can be found by adding 00767 * GeodesicExact::EllipsoidArea()/2 to the sum of \e S12 for each side of 00768 * the polygon. 00769 **********************************************************************/ 00770 Math::real EllipsoidArea() const 00771 { return 4 * Math::pi() * _c2; } 00772 ///@} 00773 00774 /** 00775 * A global instantiation of GeodesicExact with the parameters for the WGS84 00776 * ellipsoid. 00777 **********************************************************************/ 00778 static const GeodesicExact& WGS84(); 00779 00780 }; 00781 00782 } // namespace GeographicLib 00783 00784 #endif // GEOGRAPHICLIB_GEODESICEXACT_HPP