00001 /** 00002 * \file Rhumb.hpp 00003 * \brief Header for GeographicLib::Rhumb and GeographicLib::RhumbLine classes 00004 * 00005 * Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under 00006 * the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_RHUMB_HPP) 00011 #define GEOGRAPHICLIB_RHUMB_HPP 1 00012 00013 #include <GeographicLib/Constants.hpp> 00014 #include <GeographicLib/Ellipsoid.hpp> 00015 00016 namespace GeographicLib { 00017 00018 class RhumbLine; 00019 00020 /** 00021 * \brief Solve of the direct and inverse rhumb problems. 00022 * 00023 * The path of constant azimuth between two points on a ellipsoid at (\e 00024 * lat1, \e lon1) and (\e lat2, \e lon2) is called the rhumb line (also 00025 * called the loxodrome). Its length is \e s12 and its azimuth is \e azi12 00026 * and \e azi2. (The azimuth is the heading measured clockwise from north.) 00027 * 00028 * Given \e lat1, \e lon1, \e azi12, and \e s12, we can determine \e lat2, 00029 * and \e lon2. This is the \e direct rhumb problem and its solution is 00030 * given by the function Rhumb::Direct. 00031 * 00032 * Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi12 00033 * and \e s12. This is the \e inverse rhumb problem, whose solution is 00034 * given by Rhumb::Inverse. This finds the shortest such rhumb line, i.e., 00035 * the one that wraps no more than half way around the earth . 00036 * 00037 * Note that rhumb lines may be appreciably longer (up to 50%) than the 00038 * corresponding Geodesic. For example the distance between London Heathrow 00039 * and Tokyo Narita via the rhumb line is 11400 km which is 18% longer than 00040 * the geodesic distance 9600 km. 00041 * 00042 * For more information on rhumb lines see \ref rhumb. 00043 * 00044 * Example of use: 00045 * \include example-Rhumb.cpp 00046 **********************************************************************/ 00047 00048 class GEOGRAPHICLIB_EXPORT Rhumb { 00049 private: 00050 typedef Math::real real; 00051 friend class RhumbLine; 00052 Ellipsoid _ell; 00053 bool _exact; 00054 static const int tm_maxord = GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER; 00055 static inline real overflow() { 00056 // Overflow value s.t. atan(overflow_) = pi/2 00057 static const real 00058 overflow = 1 / Math::sq(std::numeric_limits<real>::epsilon()); 00059 return overflow; 00060 } 00061 static inline real tano(real x) { 00062 using std::abs; using std::tan; 00063 return 00064 2 * abs(x) == Math::pi() ? (x < 0 ? - overflow() : overflow()) : 00065 tan(x); 00066 } 00067 static inline real gd(real x) 00068 { using std::atan; using std::sinh; return atan(sinh(x)); } 00069 00070 // Use divided differences to determine (mu2 - mu1) / (psi2 - psi1) 00071 // accurately 00072 // 00073 // Definition: Df(x,y,d) = (f(x) - f(y)) / (x - y) 00074 // See: 00075 // W. M. Kahan and R. J. Fateman, 00076 // Symbolic computation of divided differences, 00077 // SIGSAM Bull. 33(3), 7-28 (1999) 00078 // http://dx.doi.org/10.1145/334714.334716 00079 // http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf 00080 00081 static inline real Dtan(real x, real y) { 00082 real d = x - y, tx = tano(x), ty = tano(y), txy = tx * ty; 00083 return d ? (2 * txy > -1 ? (1 + txy) * tano(d) : tx - ty) / d : 00084 1 + txy; 00085 } 00086 static inline real Datan(real x, real y) { 00087 using std::atan; 00088 real d = x - y, xy = x * y; 00089 return d ? (2 * xy > -1 ? atan( d / (1 + xy) ) : atan(x) - atan(y)) / d : 00090 1 / (1 + xy); 00091 } 00092 static inline real Dsin(real x, real y) { 00093 using std::sin; using std::cos; 00094 real d = (x - y)/2; 00095 return cos((x + y)/2) * (d ? sin(d) / d : 1); 00096 } 00097 static inline real Dsinh(real x, real y) { 00098 using std::sinh; using std::cosh; 00099 real d = (x - y)/2; 00100 return cosh((x + y)/2) * (d ? sinh(d) / d : 1); 00101 } 00102 static inline real Dasinh(real x, real y) { 00103 real d = x - y, 00104 hx = Math::hypot(real(1), x), hy = Math::hypot(real(1), y); 00105 return d ? Math::asinh(x*y > 0 ? d * (x + y) / (x*hy + y*hx) : 00106 x*hy - y*hx) / d : 00107 1 / hx; 00108 } 00109 static inline real Dgd(real x, real y) { 00110 using std::sinh; 00111 return Datan(sinh(x), sinh(y)) * Dsinh(x, y); 00112 } 00113 static inline real Dgdinv(real x, real y) { 00114 return Dasinh(tano(x), tano(y)) * Dtan(x, y); 00115 } 00116 // Copied from LambertConformalConic... 00117 // e * atanh(e * x) = log( ((1 + e*x)/(1 - e*x))^(e/2) ) if f >= 0 00118 // - sqrt(-e2) * atan( sqrt(-e2) * x) if f < 0 00119 inline real eatanhe(real x) const { 00120 using std::atan; 00121 return _ell._f >= 0 ? _ell._e * Math::atanh(_ell._e * x) : 00122 - _ell._e * atan(_ell._e * x); 00123 } 00124 // Copied from LambertConformalConic... 00125 // Deatanhe(x,y) = eatanhe((x-y)/(1-e^2*x*y))/(x-y) 00126 inline real Deatanhe(real x, real y) const { 00127 real t = x - y, d = 1 - _ell._e2 * x * y; 00128 return t ? eatanhe(t / d) / t : _ell._e2 / d; 00129 } 00130 // (E(x) - E(y)) / (x - y) -- E = incomplete elliptic integral of 2nd kind 00131 real DE(real x, real y) const; 00132 // (mux - muy) / (phix - phiy) using elliptic integrals 00133 real DRectifying(real latx, real laty) const; 00134 // (psix - psiy) / (phix - phiy) 00135 real DIsometric(real latx, real laty) const; 00136 00137 // (sum(c[j]*sin(2*j*x),j=1..n) - sum(c[j]*sin(2*j*x),j=1..n)) / (x - y) 00138 static real SinSeries(real x, real y, const real c[], int n); 00139 // (mux - muy) / (chix - chiy) using Krueger's series 00140 real DConformalToRectifying(real chix, real chiy) const; 00141 // (chix - chiy) / (mux - muy) using Krueger's series 00142 real DRectifyingToConformal(real mux, real muy) const; 00143 00144 // (mux - muy) / (psix - psiy) 00145 real DIsometricToRectifying(real psix, real psiy) const; 00146 // (psix - psiy) / (mux - muy) 00147 real DRectifyingToIsometric(real mux, real muy) const; 00148 00149 public: 00150 00151 /** 00152 * Constructor for a ellipsoid with 00153 * 00154 * @param[in] a equatorial radius (meters). 00155 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00156 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set 00157 * flattening to 1/\e f. 00158 * @param[in] exact if true (the default) use an addition theorem for 00159 * elliptic integrals to compute divided differences; otherwise use 00160 * series expansion (accurate for |<i>f</i>| < 0.01). 00161 * @exception GeographicErr if \e a or (1 − \e f) \e a is not 00162 * positive. 00163 * 00164 * See \ref rhumb, for a detailed description of the \e exact parameter. 00165 **********************************************************************/ 00166 Rhumb(real a, real f, bool exact = true) : _ell(a, f), _exact(exact) {} 00167 00168 /** 00169 * Solve the direct rhumb problem. 00170 * 00171 * @param[in] lat1 latitude of point 1 (degrees). 00172 * @param[in] lon1 longitude of point 1 (degrees). 00173 * @param[in] azi12 azimuth of the rhumb line (degrees). 00174 * @param[in] s12 distance between point 1 and point 2 (meters); it can be 00175 * negative. 00176 * @param[out] lat2 latitude of point 2 (degrees). 00177 * @param[out] lon2 longitude of point 2 (degrees). 00178 * 00179 * \e lat1 should be in the range [−90°, 90°]; \e lon1 and \e 00180 * azi1 should be in the range [−540°, 540°). The values of 00181 * \e lon2 and \e azi2 returned are in the range [−180°, 00182 * 180°). 00183 * 00184 * If point 1 is a pole, the cosine of its latitude is taken to be 00185 * 1/ε<sup>2</sup> (where ε is 2<sup>-52</sup>). This 00186 * position, which is extremely close to the actual pole, allows the 00187 * calculation to be carried out in finite terms. If \e s12 is large 00188 * enough that the rhumb line crosses a pole, the longitude of point 2 00189 * is indeterminate (a NaN is returned for \e lon2). 00190 **********************************************************************/ 00191 void Direct(real lat1, real lon1, real azi12, real s12, 00192 real& lat2, real& lon2) const; 00193 00194 /** 00195 * Solve the inverse rhumb problem. 00196 * 00197 * @param[in] lat1 latitude of point 1 (degrees). 00198 * @param[in] lon1 longitude of point 1 (degrees). 00199 * @param[in] lat2 latitude of point 2 (degrees). 00200 * @param[in] lon2 longitude of point 2 (degrees). 00201 * @param[out] s12 rhumb distance between point 1 and point 2 (meters). 00202 * @param[out] azi12 azimuth of the rhumb line (degrees). 00203 * 00204 * The shortest rhumb line is found. \e lat1 and \e lat2 should be in the 00205 * range [−90°, 90°]; \e lon1 and \e lon2 should be in the 00206 * range [−540°, 540°). The value of \e azi12 returned is in 00207 * the range [−180°, 180°). 00208 * 00209 * If either point is a pole, the cosine of its latitude is taken to be 00210 * 1/ε<sup>2</sup> (where ε is 2<sup>-52</sup>). This 00211 * position, which is extremely close to the actual pole, allows the 00212 * calculation to be carried out in finite terms. 00213 **********************************************************************/ 00214 void Inverse(real lat1, real lon1, real lat2, real lon2, 00215 real& s12, real& azi12) const; 00216 00217 /** 00218 * Set up to compute several points on a single rhumb line. 00219 * 00220 * @param[in] lat1 latitude of point 1 (degrees). 00221 * @param[in] lon1 longitude of point 1 (degrees). 00222 * @param[in] azi12 azimuth of the rhumb line (degrees). 00223 * @return a RhumbLine object. 00224 * 00225 * \e lat1 should be in the range [−90°, 90°]; \e lon1 and \e 00226 * azi12 should be in the range [−540°, 540°). 00227 * 00228 * If point 1 is a pole, the cosine of its latitude is taken to be 00229 * 1/ε<sup>2</sup> (where ε is 2<sup>-52</sup>). This 00230 * position, which is extremely close to the actual pole, allows the 00231 * calculation to be carried out in finite terms. 00232 **********************************************************************/ 00233 RhumbLine Line(real lat1, real lon1, real azi12) const; 00234 00235 /** \name Inspector functions. 00236 **********************************************************************/ 00237 ///@{ 00238 00239 /** 00240 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00241 * the value used in the constructor. 00242 **********************************************************************/ 00243 Math::real MajorRadius() const { return _ell.MajorRadius(); } 00244 00245 /** 00246 * @return \e f the flattening of the ellipsoid. This is the 00247 * value used in the constructor. 00248 **********************************************************************/ 00249 Math::real Flattening() const { return _ell.Flattening(); } 00250 00251 /** 00252 * A global instantiation of Rhumb with the parameters for the WGS84 00253 * ellipsoid. 00254 **********************************************************************/ 00255 static const Rhumb& WGS84(); 00256 }; 00257 00258 /** 00259 * \brief Find a sequence of points on a single rhumb line. 00260 * 00261 * RhumbLine facilitates the determination of a series of points on a single 00262 * rhumb line. The starting point (\e lat1, \e lon1) and the azimuth \e 00263 * azi12 are specified in the call to Rhumb::Line which returns a RhumbLine 00264 * object. RhumbLine.Position returns the location of point 2 a distance \e 00265 * s12 along the rhumb line. 00266 00267 * There is no public constructor for this class. (Use Rhumb::Line to create 00268 * an instance.) The Rhumb object used to create a RhumbLine must stay in 00269 * scope as long as the RhumbLine. 00270 * 00271 * Example of use: 00272 * \include example-RhumbLine.cpp 00273 **********************************************************************/ 00274 00275 class GEOGRAPHICLIB_EXPORT RhumbLine { 00276 private: 00277 typedef Math::real real; 00278 friend class Rhumb; 00279 const Rhumb& _rh; 00280 bool _exact; 00281 real _lat1, _lon1, _azi12, _salp, _calp, _mu1, _psi1, _r1; 00282 RhumbLine& operator=(const RhumbLine&); // copy assignment not allowed 00283 RhumbLine(const Rhumb& rh, real lat1, real lon1, real azi12, 00284 bool exact); 00285 public: 00286 /** 00287 * Compute the position of point 2 which is a distance \e s12 (meters) from 00288 * point 1. 00289 * 00290 * @param[in] s12 distance between point 1 and point 2 (meters); it can be 00291 * negative. 00292 * @param[out] lat2 latitude of point 2 (degrees). 00293 * @param[out] lon2 longitude of point 2 (degrees). 00294 * 00295 * The values of \e lon2 and \e azi2 returned are in the range 00296 * [−180°, 180°). 00297 * 00298 * If \e s12 is large enough that the rhumb line crosses a pole, the 00299 * longitude of point 2 is indeterminate (a NaN is returned for \e lon2). 00300 **********************************************************************/ 00301 void Position(real s12, real& lat2, real& lon2) const; 00302 00303 /** \name Inspector functions 00304 **********************************************************************/ 00305 ///@{ 00306 00307 /** 00308 * @return \e lat1 the latitude of point 1 (degrees). 00309 **********************************************************************/ 00310 Math::real Latitude() const { return _lat1; } 00311 00312 /** 00313 * @return \e lon1 the longitude of point 1 (degrees). 00314 **********************************************************************/ 00315 Math::real Longitude() const { return _lon1; } 00316 00317 /** 00318 * @return \e azi12 the azimuth of the rhumb line (degrees). 00319 **********************************************************************/ 00320 Math::real Azimuth() const { return _azi12; } 00321 00322 /** 00323 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00324 * the value inherited from the Rhumb object used in the constructor. 00325 **********************************************************************/ 00326 Math::real MajorRadius() const { return _rh.MajorRadius(); } 00327 00328 /** 00329 * @return \e f the flattening of the ellipsoid. This is the value 00330 * inherited from the Rhumb object used in the constructor. 00331 **********************************************************************/ 00332 Math::real Flattening() const { return _rh.Flattening(); } 00333 }; 00334 00335 } // namespace GeographicLib 00336 00337 #endif // GEOGRAPHICLIB_RHUMB_HPP