GeographicLib
1.21
|
00001 /** 00002 * \file Geodesic.hpp 00003 * \brief Header for GeographicLib::Geodesic class 00004 * 00005 * Copyright (c) Charles Karney (2009-2011) <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_GEODESIC_HPP) 00011 #define GEOGRAPHICLIB_GEODESIC_HPP \ 00012 "$Id: c1b085aadd7b8eabe0f9518b29531a38c152d495 $" 00013 00014 #include <GeographicLib/Constants.hpp> 00015 00016 #if !defined(GEOD_ORD) 00017 /** 00018 * The order of the expansions used by Geodesic. 00019 **********************************************************************/ 00020 #define GEOD_ORD \ 00021 (GEOGRAPHICLIB_PREC == 1 ? 6 : (GEOGRAPHICLIB_PREC == 0 ? 3 : 7)) 00022 #endif 00023 00024 namespace GeographicLib { 00025 00026 class GeodesicLine; 00027 00028 /** 00029 * \brief %Geodesic calculations 00030 * 00031 00032 * The shortest path between two points on a ellipsoid at (\e lat1, \e lon1) 00033 * and (\e lat2, \e lon2) is called the geodesic. Its length is \e s12 and 00034 * the geodesic from point 1 to point 2 has azimuths \e azi1 and \e azi2 at 00035 * the two end points. (The azimuth is the heading measured clockwise from 00036 * north. \e azi2 is the "forward" azimuth, i.e., the heading that takes you 00037 * beyond point 2 not back to point 1.) 00038 * 00039 * Given \e lat1, \e lon1, \e azi1, and \e s12, we can determine \e lat2, \e 00040 * lon2, and \e azi2. This is the \e direct geodesic problem and its 00041 * solution is given by the function Geodesic::Direct. (If \e s12 is 00042 * sufficiently large that the geodesic wraps more than halfway around the 00043 * earth, there will be another geodesic between the points with a smaller \e 00044 * s12.) 00045 * 00046 * Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi1, \e 00047 * azi2, and \e s12. This is the \e inverse geodesic problem, whose solution 00048 * is given by Geodesic::Inverse. Usually, the solution to the inverse 00049 * problem is unique. In cases where there are multiple solutions (all with 00050 * the same \e s12, of course), all the solutions can be easily generated 00051 * once a particular solution is provided. 00052 * 00053 * The standard way of specifying the direct problem is the specify the 00054 * distance \e s12 to the second point. However it is sometimes useful 00055 * instead to specify the the arc length \e a12 (in degrees) on the auxiliary 00056 * sphere. This is a mathematical construct used in solving the geodesic 00057 * problems. The solution of the direct problem in this form is provide by 00058 * Geodesic::ArcDirect. An arc length in excess of 180<sup>o</sup> indicates 00059 * that the geodesic is not a shortest path. In addition, the arc length 00060 * between an equatorial crossing and the next extremum of latitude for a 00061 * geodesic is 90<sup>o</sup>. 00062 * 00063 * This class can also calculate several other quantities related to 00064 * geodesics. These are: 00065 * - <i>reduced length</i>. If we fix the first point and increase \e azi1 00066 * by \e dazi1 (radians), the the second point is displaced \e m12 \e dazi1 00067 * in the direction \e azi2 + 90<sup>o</sup>. The quantity \e m12 is 00068 * called the "reduced length" and is symmetric under interchange of the 00069 * two points. On a curved surface the reduced length obeys a symmetry 00070 * relation, \e m12 + \e m21 = 0. On a flat surface, we have \e m12 = \e 00071 * s12. The ratio <i>s12</i>/\e m12 gives the azimuthal scale for an 00072 * azimuthal equidistant projection. 00073 * - <i>geodesic scale</i>. Consider a reference geodesic and a second 00074 * geodesic parallel to this one at point 1 and separated by a small 00075 * distance \e dt. The separation of the two geodesics at point 2 is \e 00076 * M12 \e dt where \e M12 is called the "geodesic scale". \e M21 is 00077 * defined similarly (with the geodesics being parallel at point 2). On a 00078 * flat surface, we have \e M12 = \e M21 = 1. The quantity 1/\e M12 gives 00079 * the scale of the Cassini-Soldner projection. 00080 * - <i>area</i>. Consider the quadrilateral bounded by the following lines: 00081 * the geodesic from point 1 to point 2, the meridian from point 2 to the 00082 * equator, the equator from \e lon2 to \e lon1, the meridian from the 00083 * equator to point 1. The area of this quadrilateral is represented by \e 00084 * S12 with a clockwise traversal of the perimeter counting as a positive 00085 * area and it can be used to compute the area of any simple geodesic 00086 * polygon. 00087 * 00088 * Overloaded versions of Geodesic::Direct, Geodesic::ArcDirect, and 00089 * Geodesic::Inverse allow these quantities to be returned. In addition 00090 * there are general functions Geodesic::GenDirect, and Geodesic::GenInverse 00091 * which allow an arbitrary set of results to be computed. The quantities \e 00092 * m12, \e M12, \e M21 which all specify the behavior of nearby geodesics 00093 * obey addition rules. Let points 1, 2, and 3 all lie on a single geodesic, 00094 * then 00095 * - \e m13 = \e m12 \e M23 + \e m23 \e M21 00096 * - \e M13 = \e M12 \e M23 - (1 - \e M12 \e M21) \e m23 / \e m12 00097 * - \e M31 = \e M32 \e M21 - (1 - \e M23 \e M32) \e m12 / \e m23 00098 * 00099 * Additional functionality is provided by the GeodesicLine class, which 00100 * allows a sequence of points along a geodesic to be computed. 00101 * 00102 * The calculations are accurate to better than 15 nm (15 nanometers). See 00103 * Sec. 9 of 00104 * <a href="http://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a> 00105 * for details. 00106 * 00107 * The algorithms are described in 00108 * - C. F. F. Karney, 00109 * <a href="http://arxiv.org/abs/1102.1215v1">Geodesics 00110 * on an ellipsoid of revolution</a>, 00111 * Feb. 2011; 00112 * preprint 00113 * <a href="http://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a>. 00114 * - C. F. F. Karney, 00115 * <a href="http://arxiv.org/abs/1109.4448">Algorithms for geodesics</a>, 00116 * Sept. 2011; 00117 * preprint 00118 * <a href="http://arxiv.org/abs/1109.4448">arxiv:1109.4448</a>. 00119 * . 00120 * For more information on geodesics see \ref geodesic. 00121 * 00122 * Example of use: 00123 * \include example-Geodesic.cpp 00124 * 00125 * <a href="Geod.1.html">Geod</a> is a command-line utility providing access 00126 * to the functionality of Geodesic and GeodesicLine. 00127 **********************************************************************/ 00128 00129 class GEOGRAPHIC_EXPORT Geodesic { 00130 private: 00131 typedef Math::real real; 00132 friend class GeodesicLine; 00133 static const int nA1_ = GEOD_ORD; 00134 static const int nC1_ = GEOD_ORD; 00135 static const int nC1p_ = GEOD_ORD; 00136 static const int nA2_ = GEOD_ORD; 00137 static const int nC2_ = GEOD_ORD; 00138 static const int nA3_ = GEOD_ORD; 00139 static const int nA3x_ = nA3_; 00140 static const int nC3_ = GEOD_ORD; 00141 static const int nC3x_ = (nC3_ * (nC3_ - 1)) / 2; 00142 static const int nC4_ = GEOD_ORD; 00143 static const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2; 00144 static const unsigned maxit_ = 50; 00145 00146 static const real tiny_; 00147 static const real tol0_; 00148 static const real tol1_; 00149 static const real tol2_; 00150 static const real xthresh_; 00151 00152 enum captype { 00153 CAP_NONE = 0U, 00154 CAP_C1 = 1U<<0, 00155 CAP_C1p = 1U<<1, 00156 CAP_C2 = 1U<<2, 00157 CAP_C3 = 1U<<3, 00158 CAP_C4 = 1U<<4, 00159 CAP_ALL = 0x1FU, 00160 OUT_ALL = 0x7F80U, 00161 }; 00162 00163 static real SinCosSeries(bool sinp, 00164 real sinx, real cosx, const real c[], int n) 00165 throw(); 00166 static inline real AngNormalize(real x) throw() { 00167 // Place angle in [-180, 180). Assumes x is in [-540, 540). 00168 return x >= 180 ? x - 360 : (x < -180 ? x + 360 : x); 00169 } 00170 static inline real AngRound(real x) throw() { 00171 // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57 00172 // for reals = 0.7 pm on the earth if x is an angle in degrees. (This 00173 // is about 1000 times more resolution than we get with angles around 90 00174 // degrees.) We use this to avoid having to deal with near singular 00175 // cases when x is non-zero but tiny (e.g., 1.0e-200). 00176 const real z = real(0.0625); // 1/16 00177 volatile real y = std::abs(x); 00178 // The compiler mustn't "simplify" z - (z - y) to y 00179 y = y < z ? z - (z - y) : y; 00180 return x < 0 ? -y : y; 00181 } 00182 static inline void SinCosNorm(real& sinx, real& cosx) throw() { 00183 real r = Math::hypot(sinx, cosx); 00184 sinx /= r; 00185 cosx /= r; 00186 } 00187 static real Astroid(real x, real y) throw(); 00188 00189 real _a, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2; 00190 real _A3x[nA3x_], _C3x[nC3x_], _C4x[nC4x_]; 00191 00192 void Lengths(real eps, real sig12, 00193 real ssig1, real csig1, real ssig2, real csig2, 00194 real cbet1, real cbet2, 00195 real& s12s, real& m12a, real& m0, 00196 bool scalep, real& M12, real& M21, 00197 real C1a[], real C2a[]) const throw(); 00198 real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2, 00199 real lam12, 00200 real& salp1, real& calp1, 00201 real& salp2, real& calp2, 00202 real C1a[], real C2a[]) const throw(); 00203 real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2, 00204 real salp1, real calp1, 00205 real& salp2, real& calp2, real& sig12, 00206 real& ssig1, real& csig1, real& ssig2, real& csig2, 00207 real& eps, real& domg12, bool diffp, real& dlam12, 00208 real C1a[], real C2a[], real C3a[]) 00209 const throw(); 00210 00211 // These are Maxima generated functions to provide series approximations to 00212 // the integrals for the ellipsoidal geodesic. 00213 static real A1m1f(real eps) throw(); 00214 static void C1f(real eps, real c[]) throw(); 00215 static void C1pf(real eps, real c[]) throw(); 00216 static real A2m1f(real eps) throw(); 00217 static void C2f(real eps, real c[]) throw(); 00218 00219 void A3coeff() throw(); 00220 real A3f(real eps) const throw(); 00221 void C3coeff() throw(); 00222 void C3f(real eps, real c[]) const throw(); 00223 void C4coeff() throw(); 00224 void C4f(real k2, real c[]) const throw(); 00225 00226 public: 00227 00228 /** 00229 * Bit masks for what calculations to do. These masks do double duty. 00230 * They signify to the GeodesicLine::GeodesicLine constructor and to 00231 * Geodesic::Line what capabilities should be included in the GeodesicLine 00232 * object. They also specify which results to return in the general 00233 * routines Geodesic::GenDirect and Geodesic::GenInverse routines. 00234 * GeodesicLine::mask is a duplication of this enum. 00235 **********************************************************************/ 00236 enum mask { 00237 /** 00238 * No capabilities, no output. 00239 * @hideinitializer 00240 **********************************************************************/ 00241 NONE = 0U, 00242 /** 00243 * Calculate latitude \e lat2. (It's not necessary to include this as a 00244 * capability to GeodesicLine because this is included by default.) 00245 * @hideinitializer 00246 **********************************************************************/ 00247 LATITUDE = 1U<<7 | CAP_NONE, 00248 /** 00249 * Calculate longitude \e lon2. 00250 * @hideinitializer 00251 **********************************************************************/ 00252 LONGITUDE = 1U<<8 | CAP_C3, 00253 /** 00254 * Calculate azimuths \e azi1 and \e azi2. (It's not necessary to 00255 * include this as a capability to GeodesicLine because this is included 00256 * by default.) 00257 * @hideinitializer 00258 **********************************************************************/ 00259 AZIMUTH = 1U<<9 | CAP_NONE, 00260 /** 00261 * Calculate distance \e s12. 00262 * @hideinitializer 00263 **********************************************************************/ 00264 DISTANCE = 1U<<10 | CAP_C1, 00265 /** 00266 * Allow distance \e s12 to be used as input in the direct geodesic 00267 * problem. 00268 * @hideinitializer 00269 **********************************************************************/ 00270 DISTANCE_IN = 1U<<11 | CAP_C1 | CAP_C1p, 00271 /** 00272 * Calculate reduced length \e m12. 00273 * @hideinitializer 00274 **********************************************************************/ 00275 REDUCEDLENGTH = 1U<<12 | CAP_C1 | CAP_C2, 00276 /** 00277 * Calculate geodesic scales \e M12 and \e M21. 00278 * @hideinitializer 00279 **********************************************************************/ 00280 GEODESICSCALE = 1U<<13 | CAP_C1 | CAP_C2, 00281 /** 00282 * Calculate area \e S12. 00283 * @hideinitializer 00284 **********************************************************************/ 00285 AREA = 1U<<14 | CAP_C4, 00286 /** 00287 * All capabilities. Calculate everything. 00288 * @hideinitializer 00289 **********************************************************************/ 00290 ALL = OUT_ALL| CAP_ALL, 00291 }; 00292 00293 /** \name Constructor 00294 **********************************************************************/ 00295 ///@{ 00296 /** 00297 * Constructor for a ellipsoid with 00298 * 00299 * @param[in] a equatorial radius (meters). 00300 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00301 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flattening 00302 * to 1/\e f. 00303 * 00304 * An exception is thrown if either of the axes of the ellipsoid is 00305 * non-positive. 00306 **********************************************************************/ 00307 Geodesic(real a, real f); 00308 ///@} 00309 00310 /** \name Direct geodesic problem specified in terms of distance. 00311 **********************************************************************/ 00312 ///@{ 00313 /** 00314 * Perform the direct geodesic calculation where the length of the geodesic 00315 * is specify in terms of distance. 00316 * 00317 * @param[in] lat1 latitude of point 1 (degrees). 00318 * @param[in] lon1 longitude of point 1 (degrees). 00319 * @param[in] azi1 azimuth at point 1 (degrees). 00320 * @param[in] s12 distance between point 1 and point 2 (meters); it can be 00321 * signed. 00322 * @param[out] lat2 latitude of point 2 (degrees). 00323 * @param[out] lon2 longitude of point 2 (degrees). 00324 * @param[out] azi2 (forward) azimuth at point 2 (degrees). 00325 * @param[out] m12 reduced length of geodesic (meters). 00326 * @param[out] M12 geodesic scale of point 2 relative to point 1 00327 * (dimensionless). 00328 * @param[out] M21 geodesic scale of point 1 relative to point 2 00329 * (dimensionless). 00330 * @param[out] S12 area under the geodesic (meters<sup>2</sup>). 00331 * @return \e a12 arc length of between point 1 and point 2 (degrees). 00332 * 00333 * \e lat1 should be in the range [-90, 90]; \e lon1 and \e azi1 should be 00334 * in the range [-180, 360]. The values of \e lon2 and \e azi2 returned 00335 * are in the range [-180, 180). 00336 * 00337 * If either point is at a pole, the azimuth is defined by keeping the 00338 * longitude fixed and writing \e lat = 90 - \e eps or -90 + \e eps and 00339 * taking the limit \e eps -> 0 from above. An arc length greater that 180 00340 * degrees signifies a geodesic which is not a shortest path. (For a 00341 * prolate ellipsoid, an additional condition is necessary for a shortest 00342 * path: the longitudinal extent must not exceed of 180 degrees.) 00343 * 00344 * The following functions are overloaded versions of Geodesic::Direct 00345 * which omit some of the output parameters. Note, however, that the arc 00346 * length is always computed and returned as the function value. 00347 **********************************************************************/ 00348 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00349 real& lat2, real& lon2, real& azi2, 00350 real& m12, real& M12, real& M21, real& S12) 00351 const throw() { 00352 real t; 00353 return GenDirect(lat1, lon1, azi1, false, s12, 00354 LATITUDE | LONGITUDE | AZIMUTH | 00355 REDUCEDLENGTH | GEODESICSCALE | AREA, 00356 lat2, lon2, azi2, t, m12, M12, M21, S12); 00357 } 00358 00359 /** 00360 * See the documentation for Geodesic::Direct. 00361 **********************************************************************/ 00362 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00363 real& lat2, real& lon2) 00364 const throw() { 00365 real t; 00366 return GenDirect(lat1, lon1, azi1, false, s12, 00367 LATITUDE | LONGITUDE, 00368 lat2, lon2, t, t, t, t, t, t); 00369 } 00370 00371 /** 00372 * See the documentation for Geodesic::Direct. 00373 **********************************************************************/ 00374 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00375 real& lat2, real& lon2, real& azi2) 00376 const throw() { 00377 real t; 00378 return GenDirect(lat1, lon1, azi1, false, s12, 00379 LATITUDE | LONGITUDE | AZIMUTH, 00380 lat2, lon2, azi2, t, t, t, t, t); 00381 } 00382 00383 /** 00384 * See the documentation for Geodesic::Direct. 00385 **********************************************************************/ 00386 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00387 real& lat2, real& lon2, real& azi2, real& m12) 00388 const throw() { 00389 real t; 00390 return GenDirect(lat1, lon1, azi1, false, s12, 00391 LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH, 00392 lat2, lon2, azi2, t, m12, t, t, t); 00393 } 00394 00395 /** 00396 * See the documentation for Geodesic::Direct. 00397 **********************************************************************/ 00398 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00399 real& lat2, real& lon2, real& azi2, 00400 real& M12, real& M21) 00401 const throw() { 00402 real t; 00403 return GenDirect(lat1, lon1, azi1, false, s12, 00404 LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE, 00405 lat2, lon2, azi2, t, t, M12, M21, t); 00406 } 00407 00408 /** 00409 * See the documentation for Geodesic::Direct. 00410 **********************************************************************/ 00411 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00412 real& lat2, real& lon2, real& azi2, 00413 real& m12, real& M12, real& M21) 00414 const throw() { 00415 real t; 00416 return GenDirect(lat1, lon1, azi1, false, s12, 00417 LATITUDE | LONGITUDE | AZIMUTH | 00418 REDUCEDLENGTH | GEODESICSCALE, 00419 lat2, lon2, azi2, t, m12, M12, M21, t); 00420 } 00421 ///@} 00422 00423 /** \name Direct geodesic problem specified in terms of arc length. 00424 **********************************************************************/ 00425 ///@{ 00426 /** 00427 * Perform the direct geodesic calculation where the length of the geodesic 00428 * is specify in terms of arc length. 00429 * 00430 * @param[in] lat1 latitude of point 1 (degrees). 00431 * @param[in] lon1 longitude of point 1 (degrees). 00432 * @param[in] azi1 azimuth at point 1 (degrees). 00433 * @param[in] a12 arc length between point 1 and point 2 (degrees); it can 00434 * be signed. 00435 * @param[out] lat2 latitude of point 2 (degrees). 00436 * @param[out] lon2 longitude of point 2 (degrees). 00437 * @param[out] azi2 (forward) azimuth at point 2 (degrees). 00438 * @param[out] s12 distance between point 1 and point 2 (meters). 00439 * @param[out] m12 reduced length of geodesic (meters). 00440 * @param[out] M12 geodesic scale of point 2 relative to point 1 00441 * (dimensionless). 00442 * @param[out] M21 geodesic scale of point 1 relative to point 2 00443 * (dimensionless). 00444 * @param[out] S12 area under the geodesic (meters<sup>2</sup>). 00445 * 00446 * \e lat1 should be in the range [-90, 90]; \e lon1 and \e azi1 should be 00447 * in the range [-180, 360]. The values of \e lon2 and \e azi2 returned 00448 * are in the range [-180, 180). 00449 * 00450 * If either point is at a pole, the azimuth is defined by keeping the 00451 * longitude fixed and writing \e lat = 90 - \e eps or -90 + \e eps and 00452 * taking the limit \e eps -> 0 from above. An arc length greater that 180 00453 * degrees signifies a geodesic which is not a shortest path. (For a 00454 * prolate ellipsoid, an additional condition is necessary for a shortest 00455 * path: the longitudinal extent must not exceed of 180 degrees.) 00456 * 00457 * The following functions are overloaded versions of Geodesic::Direct 00458 * which omit some of the output parameters. 00459 **********************************************************************/ 00460 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00461 real& lat2, real& lon2, real& azi2, real& s12, 00462 real& m12, real& M12, real& M21, real& S12) 00463 const throw() { 00464 GenDirect(lat1, lon1, azi1, true, a12, 00465 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE | 00466 REDUCEDLENGTH | GEODESICSCALE | AREA, 00467 lat2, lon2, azi2, s12, m12, M12, M21, S12); 00468 } 00469 00470 /** 00471 * See the documentation for Geodesic::ArcDirect. 00472 **********************************************************************/ 00473 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00474 real& lat2, real& lon2) const throw() { 00475 real t; 00476 GenDirect(lat1, lon1, azi1, true, a12, 00477 LATITUDE | LONGITUDE, 00478 lat2, lon2, t, t, t, t, t, t); 00479 } 00480 00481 /** 00482 * See the documentation for Geodesic::ArcDirect. 00483 **********************************************************************/ 00484 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00485 real& lat2, real& lon2, real& azi2) const throw() { 00486 real t; 00487 GenDirect(lat1, lon1, azi1, true, a12, 00488 LATITUDE | LONGITUDE | AZIMUTH, 00489 lat2, lon2, azi2, t, t, t, t, t); 00490 } 00491 00492 /** 00493 * See the documentation for Geodesic::ArcDirect. 00494 **********************************************************************/ 00495 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00496 real& lat2, real& lon2, real& azi2, real& s12) 00497 const throw() { 00498 real t; 00499 GenDirect(lat1, lon1, azi1, true, a12, 00500 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE, 00501 lat2, lon2, azi2, s12, t, t, t, t); 00502 } 00503 00504 /** 00505 * See the documentation for Geodesic::ArcDirect. 00506 **********************************************************************/ 00507 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00508 real& lat2, real& lon2, real& azi2, 00509 real& s12, real& m12) const throw() { 00510 real t; 00511 GenDirect(lat1, lon1, azi1, true, a12, 00512 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE | 00513 REDUCEDLENGTH, 00514 lat2, lon2, azi2, s12, m12, t, t, t); 00515 } 00516 00517 /** 00518 * See the documentation for Geodesic::ArcDirect. 00519 **********************************************************************/ 00520 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00521 real& lat2, real& lon2, real& azi2, real& s12, 00522 real& M12, real& M21) const throw() { 00523 real t; 00524 GenDirect(lat1, lon1, azi1, true, a12, 00525 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE | 00526 GEODESICSCALE, 00527 lat2, lon2, azi2, s12, t, M12, M21, t); 00528 } 00529 00530 /** 00531 * See the documentation for Geodesic::ArcDirect. 00532 **********************************************************************/ 00533 void ArcDirect(real lat1, real lon1, real azi1, real a12, 00534 real& lat2, real& lon2, real& azi2, real& s12, 00535 real& m12, real& M12, real& M21) const throw() { 00536 real t; 00537 GenDirect(lat1, lon1, azi1, true, a12, 00538 LATITUDE | LONGITUDE | AZIMUTH | DISTANCE | 00539 REDUCEDLENGTH | GEODESICSCALE, 00540 lat2, lon2, azi2, s12, m12, M12, M21, t); 00541 } 00542 ///@} 00543 00544 /** \name General version of the direct geodesic solution. 00545 **********************************************************************/ 00546 ///@{ 00547 00548 /** 00549 * The general direct geodesic calculation. Geodesic::Direct and 00550 * Geodesic::ArcDirect are defined in terms of this function. 00551 * 00552 * @param[in] lat1 latitude of point 1 (degrees). 00553 * @param[in] lon1 longitude of point 1 (degrees). 00554 * @param[in] azi1 azimuth at point 1 (degrees). 00555 * @param[in] arcmode boolean flag determining the meaning of the second 00556 * parameter. 00557 * @param[in] s12_a12 if \e arcmode is false, this is the distance between 00558 * point 1 and point 2 (meters); otherwise it is the arc length between 00559 * point 1 and point 2 (degrees); it can be signed. 00560 * @param[in] outmask a bitor'ed combination of Geodesic::mask values 00561 * specifying which of the following parameters should be set. 00562 * @param[out] lat2 latitude of point 2 (degrees). 00563 * @param[out] lon2 longitude of point 2 (degrees). 00564 * @param[out] azi2 (forward) azimuth at point 2 (degrees). 00565 * @param[out] s12 distance between point 1 and point 2 (meters). 00566 * @param[out] m12 reduced length of geodesic (meters). 00567 * @param[out] M12 geodesic scale of point 2 relative to point 1 00568 * (dimensionless). 00569 * @param[out] M21 geodesic scale of point 1 relative to point 2 00570 * (dimensionless). 00571 * @param[out] S12 area under the geodesic (meters<sup>2</sup>). 00572 * @return \e a12 arc length of between point 1 and point 2 (degrees). 00573 * 00574 * The Geodesic::mask values possible for \e outmask are 00575 * - \e outmask |= Geodesic::LATITUDE for the latitude \e lat2. 00576 * - \e outmask |= Geodesic::LONGITUDE for the latitude \e lon2. 00577 * - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2. 00578 * - \e outmask |= Geodesic::DISTANCE for the distance \e s12. 00579 * - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e 00580 * m12. 00581 * - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e 00582 * M12 and \e M21. 00583 * - \e outmask |= Geodesic::AREA for the area \e S12. 00584 * . 00585 * The function value \e a12 is always computed and returned and this 00586 * equals \e s12_a12 is \e arcmode is true. If \e outmask includes 00587 * Geodesic::DISTANCE and \e arcmode is false, then \e s12 = \e s12_a12. 00588 * It is not necessary to include Geodesic::DISTANCE_IN in \e outmask; this 00589 * is automatically included is \e arcmode is false. 00590 **********************************************************************/ 00591 Math::real GenDirect(real lat1, real lon1, real azi1, 00592 bool arcmode, real s12_a12, unsigned outmask, 00593 real& lat2, real& lon2, real& azi2, 00594 real& s12, real& m12, real& M12, real& M21, 00595 real& S12) const throw(); 00596 ///@} 00597 00598 /** \name Inverse geodesic problem. 00599 **********************************************************************/ 00600 ///@{ 00601 /** 00602 * Perform the inverse geodesic calculation. 00603 * 00604 * @param[in] lat1 latitude of point 1 (degrees). 00605 * @param[in] lon1 longitude of point 1 (degrees). 00606 * @param[in] lat2 latitude of point 2 (degrees). 00607 * @param[in] lon2 longitude of point 2 (degrees). 00608 * @param[out] s12 distance between point 1 and point 2 (meters). 00609 * @param[out] azi1 azimuth at point 1 (degrees). 00610 * @param[out] azi2 (forward) azimuth at point 2 (degrees). 00611 * @param[out] m12 reduced length of geodesic (meters). 00612 * @param[out] M12 geodesic scale of point 2 relative to point 1 00613 * (dimensionless). 00614 * @param[out] M21 geodesic scale of point 1 relative to point 2 00615 * (dimensionless). 00616 * @param[out] S12 area under the geodesic (meters<sup>2</sup>). 00617 * @return \e a12 arc length of between point 1 and point 2 (degrees). 00618 * 00619 * \e lat1 and \e lat2 should be in the range [-90, 90]; \e lon1 and \e 00620 * lon2 should be in the range [-180, 360]. The values of \e azi1 and \e 00621 * azi2 returned are in the range [-180, 180). 00622 * 00623 * If either point is at a pole, the azimuth is defined by keeping the 00624 * longitude fixed and writing \e lat = 90 - \e eps or -90 + \e eps and 00625 * taking the limit \e eps -> 0 from above. If the routine fails to 00626 * converge, then all the requested outputs are set to Math::NaN(). (Test 00627 * for such results with Math::isnan.) This is not expected to happen with 00628 * ellipsoidal models of the earth; please report all cases where this 00629 * occurs. 00630 * 00631 * The following functions are overloaded versions of Geodesic::Inverse 00632 * which omit some of the output parameters. Note, however, that the arc 00633 * length is always computed and returned as the function value. 00634 **********************************************************************/ 00635 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00636 real& s12, real& azi1, real& azi2, real& m12, 00637 real& M12, real& M21, real& S12) const throw() { 00638 return GenInverse(lat1, lon1, lat2, lon2, 00639 DISTANCE | AZIMUTH | 00640 REDUCEDLENGTH | GEODESICSCALE | AREA, 00641 s12, azi1, azi2, m12, M12, M21, S12); 00642 } 00643 00644 /** 00645 * See the documentation for Geodesic::Inverse. 00646 **********************************************************************/ 00647 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00648 real& s12) const throw() { 00649 real t; 00650 return GenInverse(lat1, lon1, lat2, lon2, 00651 DISTANCE, 00652 s12, t, t, t, t, t, t); 00653 } 00654 00655 /** 00656 * See the documentation for Geodesic::Inverse. 00657 **********************************************************************/ 00658 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00659 real& azi1, real& azi2) const throw() { 00660 real t; 00661 return GenInverse(lat1, lon1, lat2, lon2, 00662 AZIMUTH, 00663 t, azi1, azi2, t, t, t, t); 00664 } 00665 00666 /** 00667 * See the documentation for Geodesic::Inverse. 00668 **********************************************************************/ 00669 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00670 real& s12, real& azi1, real& azi2) 00671 const throw() { 00672 real t; 00673 return GenInverse(lat1, lon1, lat2, lon2, 00674 DISTANCE | AZIMUTH, 00675 s12, azi1, azi2, t, t, t, t); 00676 } 00677 00678 /** 00679 * See the documentation for Geodesic::Inverse. 00680 **********************************************************************/ 00681 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00682 real& s12, real& azi1, real& azi2, real& m12) 00683 const throw() { 00684 real t; 00685 return GenInverse(lat1, lon1, lat2, lon2, 00686 DISTANCE | AZIMUTH | REDUCEDLENGTH, 00687 s12, azi1, azi2, m12, t, t, t); 00688 } 00689 00690 /** 00691 * See the documentation for Geodesic::Inverse. 00692 **********************************************************************/ 00693 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00694 real& s12, real& azi1, real& azi2, 00695 real& M12, real& M21) const throw() { 00696 real t; 00697 return GenInverse(lat1, lon1, lat2, lon2, 00698 DISTANCE | AZIMUTH | GEODESICSCALE, 00699 s12, azi1, azi2, t, M12, M21, t); 00700 } 00701 00702 /** 00703 * See the documentation for Geodesic::Inverse. 00704 **********************************************************************/ 00705 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00706 real& s12, real& azi1, real& azi2, real& m12, 00707 real& M12, real& M21) const throw() { 00708 real t; 00709 return GenInverse(lat1, lon1, lat2, lon2, 00710 DISTANCE | AZIMUTH | 00711 REDUCEDLENGTH | GEODESICSCALE, 00712 s12, azi1, azi2, m12, M12, M21, t); 00713 } 00714 ///@} 00715 00716 /** \name General version of inverse geodesic solution. 00717 **********************************************************************/ 00718 ///@{ 00719 /** 00720 * The general inverse geodesic calculation. Geodesic::Inverse is defined 00721 * in terms of this function. 00722 * 00723 * @param[in] lat1 latitude of point 1 (degrees). 00724 * @param[in] lon1 longitude of point 1 (degrees). 00725 * @param[in] lat2 latitude of point 2 (degrees). 00726 * @param[in] lon2 longitude of point 2 (degrees). 00727 * @param[in] outmask a bitor'ed combination of Geodesic::mask values 00728 * specifying which of the following parameters should be set. 00729 * @param[out] s12 distance between point 1 and point 2 (meters). 00730 * @param[out] azi1 azimuth at point 1 (degrees). 00731 * @param[out] azi2 (forward) azimuth at point 2 (degrees). 00732 * @param[out] m12 reduced length of geodesic (meters). 00733 * @param[out] M12 geodesic scale of point 2 relative to point 1 00734 * (dimensionless). 00735 * @param[out] M21 geodesic scale of point 1 relative to point 2 00736 * (dimensionless). 00737 * @param[out] S12 area under the geodesic (meters<sup>2</sup>). 00738 * @return \e a12 arc length of between point 1 and point 2 (degrees). 00739 * 00740 * The Geodesic::mask values possible for \e outmask are 00741 * - \e outmask |= Geodesic::DISTANCE for the distance \e s12. 00742 * - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2. 00743 * - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e 00744 * m12. 00745 * - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e 00746 * M12 and \e M21. 00747 * - \e outmask |= Geodesic::AREA for the area \e S12. 00748 * . 00749 * The arc length is always computed and returned as the function value. 00750 **********************************************************************/ 00751 Math::real GenInverse(real lat1, real lon1, real lat2, real lon2, 00752 unsigned outmask, 00753 real& s12, real& azi1, real& azi2, 00754 real& m12, real& M12, real& M21, real& S12) 00755 const throw(); 00756 ///@} 00757 00758 /** \name Interface to GeodesicLine. 00759 **********************************************************************/ 00760 ///@{ 00761 00762 /** 00763 * Set up to compute several points on a singe geodesic. 00764 * 00765 * @param[in] lat1 latitude of point 1 (degrees). 00766 * @param[in] lon1 longitude of point 1 (degrees). 00767 * @param[in] azi1 azimuth at point 1 (degrees). 00768 * @param[in] caps bitor'ed combination of Geodesic::mask values 00769 * specifying the capabilities the GeodesicLine object should possess, 00770 * i.e., which quantities can be returned in calls to 00771 * GeodesicLib::Position. 00772 * 00773 * \e lat1 should be in the range [-90, 90]; \e lon1 and \e azi1 should be 00774 * in the range [-180, 360]. 00775 * 00776 * The Geodesic::mask values are 00777 * - \e caps |= Geodesic::LATITUDE for the latitude \e lat2; this is 00778 * added automatically 00779 * - \e caps |= Geodesic::LONGITUDE for the latitude \e lon2 00780 * - \e caps |= Geodesic::AZIMUTH for the latitude \e azi2; this is 00781 * added automatically 00782 * - \e caps |= Geodesic::DISTANCE for the distance \e s12 00783 * - \e caps |= Geodesic::REDUCEDLENGTH for the reduced length \e m12 00784 * - \e caps |= Geodesic::GEODESICSCALE for the geodesic scales \e M12 00785 * and \e M21 00786 * - \e caps |= Geodesic::AREA for the area \e S12 00787 * - \e caps |= Geodesic::DISTANCE_IN permits the length of the 00788 * geodesic to be given in terms of \e s12; without this capability the 00789 * length can only be specified in terms of arc length. 00790 * . 00791 * The default value of \e caps is Geodesic::ALL which turns on all the 00792 * capabilities. 00793 * 00794 * If the point is at a pole, the azimuth is defined by keeping the \e lon1 00795 * fixed and writing \e lat1 = 90 - \e eps or -90 + \e eps and taking the 00796 * limit \e eps -> 0 from above. 00797 **********************************************************************/ 00798 GeodesicLine Line(real lat1, real lon1, real azi1, unsigned caps = ALL) 00799 const throw(); 00800 00801 ///@} 00802 00803 /** \name Inspector functions. 00804 **********************************************************************/ 00805 ///@{ 00806 00807 /** 00808 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00809 * the value used in the constructor. 00810 **********************************************************************/ 00811 Math::real MajorRadius() const throw() { return _a; } 00812 00813 /** 00814 * @return \e f the flattening of the ellipsoid. This is the 00815 * value used in the constructor. 00816 **********************************************************************/ 00817 Math::real Flattening() const throw() { return _f; } 00818 00819 /// \cond SKIP 00820 /** 00821 * <b>DEPRECATED</b> 00822 * @return \e r the inverse flattening of the ellipsoid. 00823 **********************************************************************/ 00824 Math::real InverseFlattening() const throw() { return 1/_f; } 00825 /// \endcond 00826 00827 /** 00828 * @return total area of ellipsoid in meters<sup>2</sup>. The area of a 00829 * polygon encircling a pole can be found by adding 00830 * Geodesic::EllipsoidArea()/2 to the sum of \e S12 for each side of the 00831 * polygon. 00832 **********************************************************************/ 00833 Math::real EllipsoidArea() const throw() 00834 { return 4 * Math::pi<real>() * _c2; } 00835 ///@} 00836 00837 /** 00838 * A global instantiation of Geodesic with the parameters for the WGS84 00839 * ellipsoid. 00840 **********************************************************************/ 00841 static const Geodesic WGS84; 00842 00843 }; 00844 00845 } // namespace GeographicLib 00846 00847 #endif // GEOGRAPHICLIB_GEODESIC_HPP