00001 /** 00002 * \file Geodesic.hpp 00003 * \brief Header for GeographicLib::Geodesic class 00004 * 00005 * Copyright (c) Charles Karney (2009, 2010) <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 6827 2010-05-20 19:56:18Z 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 * The shortest path between two points on a ellipsoid at (\e lat1, \e lon1) 00030 * and (\e lat2, \e lon2) is called the geodesic. Its length is \e s12 and 00031 * the geodesic from point 1 to point 2 has azimuths \e azi1 and \e azi2 at 00032 * the two end points. (The azimuth is the heading measured clockwise from 00033 * north. \e azi2 is the "forward" azimuth, i.e., the heading that takes you 00034 * beyond point 2 not back to point 1.) 00035 * 00036 * If we fix the first point and increase \e s12 by \e ds12, then the second 00037 * point is displaced \e ds12 in the direction \e azi2. Similarly we 00038 * increase \e azi1 by \e dazi1 (radians), the the second point is displaced 00039 * \e m12 \e dazi1 in the direction \e azi2 + 90<sup>o</sup>. The quantity 00040 * \e m12 is called the "reduced length" and is symmetric under interchange 00041 * of the two points. On a flat surface, he have \e m12 = \e s12. The ratio 00042 * \e s12/\e m12 gives the azimuthal scale for an azimuthal equidistant 00043 * projection. 00044 * 00045 * Given \e lat1, \e lon1, \e azi1, and \e s12, we can determine \e lat2, \e 00046 * lon2, \e azi2, \e m12. This is the \e direct geodesic problem. (If \e 00047 * s12 is sufficiently large that the geodesic wraps more than halfway around 00048 * the earth, there will be another geodesic between the points with a 00049 * smaller \e s12.) 00050 * 00051 * Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi1, \e 00052 * azi2, \e s12, \e m12. This is the \e inverse geodesic problem. Usually, 00053 * the solution to the inverse problem is unique. In cases where there are 00054 * muliple solutions (all with the same \e s12, of course), all the solutions 00055 * can be easily generated once a particular solution is provided. 00056 * 00057 * As an alternative to using distance to measure \e s12, the class can also 00058 * use the arc length \e a12 (in degrees) on the auxiliary sphere. This is a 00059 * mathematical construct used in solving the geodesic problems. However, an 00060 * arc length in excess of 180<sup>o</sup> indicates that the geodesic is not 00061 * a shortest path. In addition, the arc length between an equatorial 00062 * crossing and the next extremum of latitude for a geodesic is 00063 * 90<sup>o</sup>. 00064 * 00065 * The calculations are accurate to better than 15 nm. (See \ref geoderrors 00066 * for details.) 00067 * 00068 * For more information on geodesics see \ref geodesic. 00069 **********************************************************************/ 00070 00071 class Geodesic { 00072 private: 00073 typedef Math::real real; 00074 friend class GeodesicLine; 00075 static const int nA1 = GEOD_ORD, nC1 = GEOD_ORD, nC1p = GEOD_ORD, 00076 nA2 = GEOD_ORD, nC2 = GEOD_ORD, nA3 = GEOD_ORD, nC3 = GEOD_ORD; 00077 static const unsigned maxit = 50; 00078 00079 static inline real sq(real x) throw() { return x * x; } 00080 void Lengths(real eps, real sig12, 00081 real ssig1, real csig1, real ssig2, real csig2, 00082 real cbet1, real cbet2, 00083 real& s12s, real& m12a, real& m0, 00084 real tc[], real zc[]) const throw(); 00085 static real Astroid(real R, real z) throw(); 00086 real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2, 00087 real lam12, 00088 real& salp1, real& calp1, 00089 real& salp2, real& calp2, 00090 real C1a[], real C2a[]) const throw(); 00091 real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2, 00092 real salp1, real calp1, 00093 real& salp2, real& calp2, real& sig12, 00094 real& ssig1, real& csig1, real& ssig2, real& csig2, 00095 real& eps, bool diffp, real& dlam12, 00096 real C1a[], real C2a[], real C3a[]) 00097 const throw(); 00098 00099 static const real eps2, tol0, tol1, tol2, xthresh; 00100 const real _a, _r, _f, _f1, _e2, _ep2, _n, _b, _etol2; 00101 static real SinSeries(real sinx, real cosx, const real c[], int n) 00102 throw(); 00103 static inline real AngNormalize(real x) throw() { 00104 // Place angle in [-180, 180). Assumes x is in [-540, 540). 00105 return x >= 180 ? x - 360 : x < -180 ? x + 360 : x; 00106 } 00107 static inline real AngRound(real x) throw() { 00108 // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57 00109 // for reals = 0.7 pm on the earth if x is an angle in degrees. (This 00110 // is about 1000 times more resolution than we get with angles around 90 00111 // degrees.) We use this to avoid having to deal with near singular 00112 // cases when x is non-zero but tiny (e.g., 1.0e-200). 00113 const real z = real(0.0625); // 1/16 00114 volatile real y = std::abs(x); 00115 // The compiler mustn't "simplify" z - (z - y) to y 00116 y = y < z ? z - (z - y) : y; 00117 return x < 0 ? -y : y; 00118 } 00119 static inline void SinCosNorm(real& sinx, real& cosx) throw() { 00120 real r = Math::hypot(sinx, cosx); 00121 sinx /= r; 00122 cosx /= r; 00123 } 00124 // These are Maxima generated functions to provide series approximations to 00125 // the integrals for the ellipsoidal geodesic. 00126 static real A1m1f(real eps) throw(); 00127 static void C1f(real eps, real c[]) throw(); 00128 static void C1pf(real eps, real c[]) throw(); 00129 static real A2m1f(real eps) throw(); 00130 static void C2f(real eps, real c[]) throw(); 00131 static real A3f(real f, real eps) throw(); 00132 static void C3f(real f, real eps, real c[]) throw(); 00133 00134 public: 00135 00136 /** 00137 * Constructor for a ellipsoid with equatorial radius \e a (meters) and 00138 * reciprocal flattening \e r. Setting \e r = 0 implies \e r = inf or 00139 * flattening = 0 (i.e., a sphere). Negative \e r indicates a prolate 00140 * ellipsoid. An exception is thrown if \e a is not positive. 00141 **********************************************************************/ 00142 Geodesic(real a, real r); 00143 00144 /** 00145 * Perform the direct geodesic calculation. Given a latitude, \e lat1, 00146 * longitude, \e lon1, and azimuth \e azi1 (degrees) for point 1 and a 00147 * range, \e s12 (meters) from point 1 to point 2, return the latitude, \e 00148 * lat2, longitude, \e lon2, and forward azimuth, \e azi2 (degrees) for 00149 * point 2 and the reduced length \e m12 (meters). If \e arcmode (default 00150 * false) is set to true, \e s12 is interpreted as the arc length \e a12 00151 * (degrees) on the auxiliary sphere. An arc length greater that 180 00152 * degrees results in a geodesic which is not a shortest path. For a 00153 * prolate ellipsoid, an additional condition is necessary for a shortest 00154 * path: the longitudinal extent must not exceed of 180 degrees. Returned 00155 * value is the arc length \e a12 (degrees) if \e arcmode is false, 00156 * otherwise it is the distance \e s12 (meters) 00157 **********************************************************************/ 00158 Math::real Direct(real lat1, real lon1, real azi1, real s12, 00159 real& lat2, real& lon2, real& azi2, real& m12, 00160 bool arcmode = false) const throw(); 00161 00162 /** 00163 * Perform the inverse geodesic calculation. Given a latitude, \e lat1, 00164 * longitude, \e lon1, for point 1 and a latitude, \e lat2, longitude, \e 00165 * lon2, for point 2 (all in degrees), return the geodesic distance, \e s12 00166 * (meters), the forward azimuths, \e azi1 and \e azi2 (degrees) at points 00167 * 1 and 2, and the reduced length \e m12 (meters). Returned value is the 00168 * arc length \e a12 (degrees) on the auxiliary sphere. The routine uses 00169 * an iterative method. If the method fails to converge, the negative of 00170 * the distances (\e s12, \e m12, and \e a12) and reverse of the azimuths 00171 * are returned. This is not expected to happen with ellipsoidal models of 00172 * the earth. Please report all cases where this occurs. 00173 **********************************************************************/ 00174 Math::real Inverse(real lat1, real lon1, real lat2, real lon2, 00175 real& s12, real& azi1, real& azi2, real& m12) 00176 const throw(); 00177 00178 /** 00179 * Set up to do a series of ranges. This returns a GeodesicLine object 00180 * with point 1 given by latitude, \e lat1, longitude, \e lon1, and azimuth 00181 * \e azi1 (degrees). Calls to GeodesicLine::Position return the 00182 * position and azimuth for point 2 a specified distance away. Using 00183 * GeodesicLine::Position is approximately 2.1 faster than calling 00184 * Geodesic::Direct. 00185 **********************************************************************/ 00186 GeodesicLine Line(real lat1, real lon1, real azi1) 00187 const throw(); 00188 00189 /** 00190 * The major radius of the ellipsoid (meters). This is that value of \e a 00191 * used in the constructor. 00192 **********************************************************************/ 00193 Math::real MajorRadius() const throw() { return _a; } 00194 00195 /** 00196 * The inverse flattening of the ellipsoid. This is that value of \e r 00197 * used in the constructor. A value of 0 is returned for a sphere 00198 * (infinite inverse flattening). 00199 **********************************************************************/ 00200 Math::real InverseFlattening() const throw() { return _r; } 00201 00202 /** 00203 * A global instantiation of Geodesic with the parameters for the WGS84 00204 * ellipsoid. 00205 **********************************************************************/ 00206 const static Geodesic WGS84; 00207 00208 }; 00209 00210 /** 00211 * \brief A geodesic line. 00212 * 00213 * GeodesicLine facilitates the determination of a series of points on a 00214 * single geodesic. Geodesic.Line returns a GeodesicLine object with the 00215 * geodesic defined by by \e lat1, \e lon1, and \e azi1. 00216 * GeodesicLine.Position returns the \e lat2, \e lon2, \e azi2, and \e m12 00217 * given \e s12. An example of use of this class is: 00218 \verbatim 00219 // Print positions on a geodesic going through latitude 30, 00220 // longitude 10 at azimuth 80. Points at intervals of 10km 00221 // in the range [-1000km, 1000km] are given. 00222 GeographicLib::GeodesicLine 00223 line(GeographicLib::Geodesic::WGS84.Line(30.0, 10.0, 80.0)); 00224 double step = 10e3; 00225 for (int s = -100; s <= 100; ++s) { 00226 double lat2, lon2, azi2, m12; 00227 double s12 = s * step; 00228 line.Position(s12, lat2, lon2, azi2, m12); 00229 cout << s12 << " " << lat2 << " " << lon2 << " " 00230 << azi2 << " " << m12 << "\n"; 00231 } 00232 \endverbatim 00233 * The default copy constructor and assignment operators work with this 00234 * class, so that, for example, the previous example could start 00235 \verbatim 00236 GeographicLib::GeodesicLine line; 00237 line = GeographicLib::Geodesic::WGS84.Line(30.0, 10.0, 80.0); 00238 ... 00239 \endverbatim 00240 * Similarly, a vector can be used to hold GeodesicLine objects. 00241 * 00242 * The calculations are accurate to better than 12 nm. (See \ref geoderrors 00243 * for details.) 00244 **********************************************************************/ 00245 00246 class GeodesicLine { 00247 private: 00248 typedef Math::real real; 00249 friend class Geodesic; 00250 static const int nC1 = Geodesic::nC1, nC1p = Geodesic::nC1p, 00251 nC2 = Geodesic::nC2, nC3 = Geodesic::nC3; 00252 00253 real _lat1, _lon1, _azi1; 00254 real _a, _r, _b, _f1, _salp0, _calp0, _k2, 00255 _ssig1, _csig1, _stau1, _ctau1, _somg1, _comg1, 00256 _A1m1, _A2m1, _A3c, _B11, _B21, _B31; 00257 // index zero elements of these arrays are unused 00258 real _C1a[nC1 + 1], _C1pa[nC1p + 1], _C2a[nC2 + 1], _C3a[nC3]; 00259 00260 static inline real sq(real x) throw() { return x * x; } 00261 GeodesicLine(const Geodesic& g, real lat1, real lon1, real azi1) 00262 throw(); 00263 public: 00264 00265 /** 00266 * A default constructor. If GeodesicLine::Position is called on the 00267 * resulting object, it returns immediately (without doing any 00268 * calculations). The object should be set with a call to Geodesic::Line. 00269 * Use Init() to test whether object is still in this uninitialized state. 00270 **********************************************************************/ 00271 GeodesicLine() throw() : _b(0) {}; 00272 00273 /** 00274 * Return the latitude, \e lat2, longitude, \e lon2, and forward azimuth, 00275 * \e azi2 (degrees) of the point 2 which is a distance, \e s12 (in 00276 * meters), from point 1. Also return the reduced length \e m12 (meters). 00277 * \e s12 can be signed. If \e arcmode (default false) is set to true, \e 00278 * s12 is interpreted as the arc length \e a12 (in degrees) on the 00279 * auxiliary sphere. Returned value is the arc length \e a12 (degrees) if 00280 * \e arcmode is false, otherwise it is the distance \e s12 (meters). 00281 **********************************************************************/ 00282 Math::real Position(real s12, real& lat2, real& lon2, 00283 real& azi2, real &m12, bool arcmode = false) 00284 const throw(); 00285 00286 /** 00287 * Return the scale of the geodesic line extending an arc length \e a12 00288 * (degrees) from point 1 to point 2. \e M12 (a number) measures the 00289 * convergence of initially parallel geodesics. It is defined by the 00290 * following construction: starting at point 1 proceed at azimuth \e azi1 + 00291 * 90<sup>o</sup> a small distance \e dt; turn -90<sup>o</sup> and proceed 00292 * a distance \e s12 (\e not the arc length \e a12); the distance to point 00293 * 2 is given by \e M12 \e dt. \e M21 is defined analogously. 00294 **********************************************************************/ 00295 void Scale(real a12, real& M12, real& M21) const throw(); 00296 00297 /** 00298 * Has this object been initialized so that Position can be called? 00299 **********************************************************************/ 00300 bool Init() const throw() { return _b > 0; } 00301 00302 /** 00303 * Return the latitude of point 1 (degrees). 00304 **********************************************************************/ 00305 Math::real Latitude() const throw() { return Init() ? _lat1 : 0; } 00306 00307 /** 00308 * Return the longitude of point 1 (degrees). 00309 **********************************************************************/ 00310 Math::real Longitude() const throw() { return Init() ? _lon1 : 0; } 00311 00312 /** 00313 * Return the azimuth (degrees) of the geodesic line as it passes through 00314 * point 1. 00315 **********************************************************************/ 00316 Math::real Azimuth() const throw() { return Init() ? _azi1 : 0; } 00317 00318 /** 00319 * Return the azimuth (degrees) of the geodesic line as it crosses the 00320 * equator in a northward direction. 00321 **********************************************************************/ 00322 Math::real EquatorialAzimuth() const throw() { 00323 return Init() ? atan2(_salp0, _calp0) / Constants::degree() : 0; 00324 } 00325 00326 /** 00327 * Return the arc length (degrees) between the northward equatorial 00328 * crossing and point 1. 00329 **********************************************************************/ 00330 Math::real EquatorialArc() const throw() { 00331 return Init() ? atan2(_ssig1, _csig1) / Constants::degree() : 0; 00332 } 00333 00334 /** 00335 * The major radius of the ellipsoid (meters). This is that value of \e a 00336 * inherited from the Geodesic object used in the constructor. 00337 **********************************************************************/ 00338 Math::real MajorRadius() const throw() { return Init() ? _a : 0; } 00339 00340 /** 00341 * The inverse flattening of the ellipsoid. This is that value of \e r 00342 * inherited from the Geodesic object used in the constructor. A value of 00343 * 0 is returned for a sphere (infinite inverse flattening). 00344 **********************************************************************/ 00345 Math::real InverseFlattening() const throw() { return Init() ? _r : 0; } 00346 }; 00347 00348 } // namespace GeographicLib 00349 #endif