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