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