/** * \file Rhumb.hpp * \brief Header for GeographicLib::Rhumb and GeographicLib::RhumbLine classes * * Copyright (c) Charles Karney (2014-2015) and licensed * under the MIT/X11 License. For more information, see * http://geographiclib.sourceforge.net/ **********************************************************************/ #if !defined(GEOGRAPHICLIB_RHUMB_HPP) #define GEOGRAPHICLIB_RHUMB_HPP 1 #include #include #if !defined(GEOGRAPHICLIB_RHUMBAREA_ORDER) /** * The order of the series approximation used in rhumb area calculations. * GEOGRAPHICLIB_RHUMBAREA_ORDER can be set to any integer in [4, 8]. **********************************************************************/ # define GEOGRAPHICLIB_RHUMBAREA_ORDER \ (GEOGRAPHICLIB_PRECISION == 2 ? 6 : \ (GEOGRAPHICLIB_PRECISION == 1 ? 4 : 8)) #endif namespace GeographicLib { class RhumbLine; template class PolygonAreaT; /** * \brief Solve of the direct and inverse rhumb problems. * * The path of constant azimuth between two points on a ellipsoid at (\e * lat1, \e lon1) and (\e lat2, \e lon2) is called the rhumb line (also * called the loxodrome). Its length is \e s12 and its azimuth is \e azi12. * (The azimuth is the heading measured clockwise from north.) * * Given \e lat1, \e lon1, \e azi12, and \e s12, we can determine \e lat2, * and \e lon2. This is the \e direct rhumb problem and its solution is * given by the function Rhumb::Direct. * * Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi12 * and \e s12. This is the \e inverse rhumb problem, whose solution is given * by Rhumb::Inverse. This finds the shortest such rhumb line, i.e., the one * that wraps no more than half way around the earth. If the end points are * on opposite meridians, there are two shortest rhumb lines and the * east-going one is chosen. * * These routines also optionally calculate the area under the rhumb line, \e * S12. This is the area, measured counter-clockwise, of the rhumb line * quadrilateral with corners (lat1,lon1), (0,lon1), * (0,lon2), and (lat2,lon2). * * Note that rhumb lines may be appreciably longer (up to 50%) than the * corresponding Geodesic. For example the distance between London Heathrow * and Tokyo Narita via the rhumb line is 11400 km which is 18% longer than * the geodesic distance 9600 km. * * For more information on rhumb lines see \ref rhumb. * * Example of use: * \include example-Rhumb.cpp **********************************************************************/ class GEOGRAPHICLIB_EXPORT Rhumb { private: typedef Math::real real; friend class RhumbLine; template friend class PolygonAreaT; Ellipsoid _ell; bool _exact; real _c2; static const int tm_maxord = GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER; static const int maxpow_ = GEOGRAPHICLIB_RHUMBAREA_ORDER; // _R[0] unused real _R[maxpow_ + 1]; static inline real gd(real x) { using std::atan; using std::sinh; return atan(sinh(x)); } // Use divided differences to determine (mu2 - mu1) / (psi2 - psi1) // accurately // // Definition: Df(x,y,d) = (f(x) - f(y)) / (x - y) // See: // W. M. Kahan and R. J. Fateman, // Symbolic computation of divided differences, // SIGSAM Bull. 33(3), 7-28 (1999) // https://dx.doi.org/10.1145/334714.334716 // http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf static inline real Dlog(real x, real y) { real t = x - y; return t ? 2 * Math::atanh(t / (x + y)) / t : 1 / x; } // N.B., x and y are in degrees static inline real Dtan(real x, real y) { real d = x - y, tx = Math::tand(x), ty = Math::tand(y), txy = tx * ty; return d ? (2 * txy > -1 ? (1 + txy) * Math::tand(d) : tx - ty) / (d * Math::degree()) : 1 + txy; } static inline real Datan(real x, real y) { using std::atan; real d = x - y, xy = x * y; return d ? (2 * xy > -1 ? atan( d / (1 + xy) ) : atan(x) - atan(y)) / d : 1 / (1 + xy); } static inline real Dsin(real x, real y) { using std::sin; using std::cos; real d = (x - y) / 2; return cos((x + y)/2) * (d ? sin(d) / d : 1); } static inline real Dsinh(real x, real y) { using std::sinh; using std::cosh; real d = (x - y) / 2; return cosh((x + y) / 2) * (d ? sinh(d) / d : 1); } static inline real Dcosh(real x, real y) { using std::sinh; real d = (x - y) / 2; return sinh((x + y) / 2) * (d ? sinh(d) / d : 1); } static inline real Dasinh(real x, real y) { real d = x - y, hx = Math::hypot(real(1), x), hy = Math::hypot(real(1), y); return d ? Math::asinh(x*y > 0 ? d * (x + y) / (x*hy + y*hx) : x*hy - y*hx) / d : 1 / hx; } static inline real Dgd(real x, real y) { using std::sinh; return Datan(sinh(x), sinh(y)) * Dsinh(x, y); } // N.B., x and y are the tangents of the angles static inline real Dgdinv(real x, real y) { return Dasinh(x, y) / Datan(x, y); } // Copied from LambertConformalConic... // Deatanhe(x,y) = eatanhe((x-y)/(1-e^2*x*y))/(x-y) inline real Deatanhe(real x, real y) const { real t = x - y, d = 1 - _ell._e2 * x * y; return t ? Math::eatanhe(t / d, _ell._es) / t : _ell._e2 / d; } // (E(x) - E(y)) / (x - y) -- E = incomplete elliptic integral of 2nd kind real DE(real x, real y) const; // (mux - muy) / (phix - phiy) using elliptic integrals real DRectifying(real latx, real laty) const; // (psix - psiy) / (phix - phiy) real DIsometric(real latx, real laty) const; // (sum(c[j]*sin(2*j*x),j=1..n) - sum(c[j]*sin(2*j*x),j=1..n)) / (x - y) static real SinCosSeries(bool sinp, real x, real y, const real c[], int n); // (mux - muy) / (chix - chiy) using Krueger's series real DConformalToRectifying(real chix, real chiy) const; // (chix - chiy) / (mux - muy) using Krueger's series real DRectifyingToConformal(real mux, real muy) const; // (mux - muy) / (psix - psiy) // N.B., psix and psiy are in degrees real DIsometricToRectifying(real psix, real psiy) const; // (psix - psiy) / (mux - muy) real DRectifyingToIsometric(real mux, real muy) const; real MeanSinXi(real psi1, real psi2) const; // The following two functions (with lots of ignored arguments) mimic the // interface to the corresponding Geodesic function. These are needed by // PolygonAreaT. void GenDirect(real lat1, real lon1, real azi12, bool, real s12, unsigned outmask, real& lat2, real& lon2, real&, real&, real&, real&, real&, real& S12) const { GenDirect(lat1, lon1, azi12, s12, outmask, lat2, lon2, S12); } void GenInverse(real lat1, real lon1, real lat2, real lon2, unsigned outmask, real& s12, real& azi12, real&, real& , real& , real& , real& S12) const { GenInverse(lat1, lon1, lat2, lon2, outmask, s12, azi12, S12); } public: /** * Bit masks for what calculations to do. They specify which results to * return in the general routines Rhumb::GenDirect and Rhumb::GenInverse * routines. RhumbLine::mask is a duplication of this enum. **********************************************************************/ enum mask { /** * No output. * @hideinitializer **********************************************************************/ NONE = 0U, /** * Calculate latitude \e lat2. * @hideinitializer **********************************************************************/ LATITUDE = 1U<<7, /** * Calculate longitude \e lon2. * @hideinitializer **********************************************************************/ LONGITUDE = 1U<<8, /** * Calculate azimuth \e azi12. * @hideinitializer **********************************************************************/ AZIMUTH = 1U<<9, /** * Calculate distance \e s12. * @hideinitializer **********************************************************************/ DISTANCE = 1U<<10, /** * Calculate area \e S12. * @hideinitializer **********************************************************************/ AREA = 1U<<14, /** * Unroll \e lon2 in the direct calculation. * @hideinitializer **********************************************************************/ LONG_UNROLL = 1U<<15, /** * Calculate everything. (LONG_UNROLL is not included in this mask.) * @hideinitializer **********************************************************************/ ALL = 0x7F80U, }; /** * Constructor for a ellipsoid with * * @param[in] a equatorial radius (meters). * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. * Negative \e f gives a prolate ellipsoid. * @param[in] exact if true (the default) use an addition theorem for * elliptic integrals to compute divided differences; otherwise use * series expansion (accurate for |f| < 0.01). * @exception GeographicErr if \e a or (1 − \e f) \e a is not * positive. * * See \ref rhumb, for a detailed description of the \e exact parameter. **********************************************************************/ Rhumb(real a, real f, bool exact = true); /** * Solve the direct rhumb problem returning also the area. * * @param[in] lat1 latitude of point 1 (degrees). * @param[in] lon1 longitude of point 1 (degrees). * @param[in] azi12 azimuth of the rhumb line (degrees). * @param[in] s12 distance between point 1 and point 2 (meters); it can be * negative. * @param[out] lat2 latitude of point 2 (degrees). * @param[out] lon2 longitude of point 2 (degrees). * @param[out] S12 area under the rhumb line (meters²). * * \e lat1 should be in the range [−90°, 90°]. The value of * \e lon2 returned is in the range [−180°, 180°). * * If point 1 is a pole, the cosine of its latitude is taken to be * 1/ε² (where ε is 2^-52). This * position, which is extremely close to the actual pole, allows the * calculation to be carried out in finite terms. If \e s12 is large * enough that the rhumb line crosses a pole, the longitude of point 2 * is indeterminate (a NaN is returned for \e lon2 and \e S12). **********************************************************************/ void Direct(real lat1, real lon1, real azi12, real s12, real& lat2, real& lon2, real& S12) const { GenDirect(lat1, lon1, azi12, s12, LATITUDE | LONGITUDE | AREA, lat2, lon2, S12); } /** * Solve the direct rhumb problem without the area. **********************************************************************/ void Direct(real lat1, real lon1, real azi12, real s12, real& lat2, real& lon2) const { real t; GenDirect(lat1, lon1, azi12, s12, LATITUDE | LONGITUDE, lat2, lon2, t); } /** * The general direct rhumb problem. Rhumb::Direct is defined in terms * of this function. * * @param[in] lat1 latitude of point 1 (degrees). * @param[in] lon1 longitude of point 1 (degrees). * @param[in] azi12 azimuth of the rhumb line (degrees). * @param[in] s12 distance between point 1 and point 2 (meters); it can be * negative. * @param[in] outmask a bitor'ed combination of Rhumb::mask values * specifying which of the following parameters should be set. * @param[out] lat2 latitude of point 2 (degrees). * @param[out] lon2 longitude of point 2 (degrees). * @param[out] S12 area under the rhumb line (meters²). * * The Rhumb::mask values possible for \e outmask are * - \e outmask |= Rhumb::LATITUDE for the latitude \e lat2; * - \e outmask |= Rhumb::LONGITUDE for the latitude \e lon2; * - \e outmask |= Rhumb::AREA for the area \e S12; * - \e outmask |= Rhumb::ALL for all of the above; * - \e outmask |= Rhumb::LONG_UNROLL to unroll \e lon2 instead of wrapping * it into the range [−180°, 180°). * . * With the Rhumb::LONG_UNROLL bit set, the quantity \e lon2 − * \e lon1 indicates how many times and in what sense the rhumb line * encircles the ellipsoid. **********************************************************************/ void GenDirect(real lat1, real lon1, real azi12, real s12, unsigned outmask, real& lat2, real& lon2, real& S12) const; /** * Solve the inverse rhumb problem returning also the area. * * @param[in] lat1 latitude of point 1 (degrees). * @param[in] lon1 longitude of point 1 (degrees). * @param[in] lat2 latitude of point 2 (degrees). * @param[in] lon2 longitude of point 2 (degrees). * @param[out] s12 rhumb distance between point 1 and point 2 (meters). * @param[out] azi12 azimuth of the rhumb line (degrees). * @param[out] S12 area under the rhumb line (meters²). * * The shortest rhumb line is found. If the end points are on opposite * meridians, there are two shortest rhumb lines and the east-going one is * chosen. \e lat1 and \e lat2 should be in the range [−90°, * 90°]. The value of \e azi12 returned is in the range * [−180°, 180°). * * If either point is a pole, the cosine of its latitude is taken to be * 1/ε² (where ε is 2^-52). This * position, which is extremely close to the actual pole, allows the * calculation to be carried out in finite terms. **********************************************************************/ void Inverse(real lat1, real lon1, real lat2, real lon2, real& s12, real& azi12, real& S12) const { GenInverse(lat1, lon1, lat2, lon2, DISTANCE | AZIMUTH | AREA, s12, azi12, S12); } /** * Solve the inverse rhumb problem without the area. **********************************************************************/ void Inverse(real lat1, real lon1, real lat2, real lon2, real& s12, real& azi12) const { real t; GenInverse(lat1, lon1, lat2, lon2, DISTANCE | AZIMUTH, s12, azi12, t); } /** * The general inverse rhumb problem. Rhumb::Inverse is defined in terms * of this function. * * @param[in] lat1 latitude of point 1 (degrees). * @param[in] lon1 longitude of point 1 (degrees). * @param[in] lat2 latitude of point 2 (degrees). * @param[in] lon2 longitude of point 2 (degrees). * @param[in] outmask a bitor'ed combination of Rhumb::mask values * specifying which of the following parameters should be set. * @param[out] s12 rhumb distance between point 1 and point 2 (meters). * @param[out] azi12 azimuth of the rhumb line (degrees). * @param[out] S12 area under the rhumb line (meters²). * * The Rhumb::mask values possible for \e outmask are * - \e outmask |= Rhumb::DISTANCE for the latitude \e s12; * - \e outmask |= Rhumb::AZIMUTH for the latitude \e azi12; * - \e outmask |= Rhumb::AREA for the area \e S12; * - \e outmask |= Rhumb::ALL for all of the above; **********************************************************************/ void GenInverse(real lat1, real lon1, real lat2, real lon2, unsigned outmask, real& s12, real& azi12, real& S12) const; /** * Set up to compute several points on a single rhumb line. * * @param[in] lat1 latitude of point 1 (degrees). * @param[in] lon1 longitude of point 1 (degrees). * @param[in] azi12 azimuth of the rhumb line (degrees). * @return a RhumbLine object. * * \e lat1 should be in the range [−90°, 90°]. * * If point 1 is a pole, the cosine of its latitude is taken to be * 1/ε² (where ε is 2^-52). This * position, which is extremely close to the actual pole, allows the * calculation to be carried out in finite terms. **********************************************************************/ RhumbLine Line(real lat1, real lon1, real azi12) const; /** \name Inspector functions. **********************************************************************/ ///@{ /** * @return \e a the equatorial radius of the ellipsoid (meters). This is * the value used in the constructor. **********************************************************************/ Math::real MajorRadius() const { return _ell.MajorRadius(); } /** * @return \e f the flattening of the ellipsoid. This is the * value used in the constructor. **********************************************************************/ Math::real Flattening() const { return _ell.Flattening(); } Math::real EllipsoidArea() const { return _ell.Area(); } /** * A global instantiation of Rhumb with the parameters for the WGS84 * ellipsoid. **********************************************************************/ static const Rhumb& WGS84(); }; /** * \brief Find a sequence of points on a single rhumb line. * * RhumbLine facilitates the determination of a series of points on a single * rhumb line. The starting point (\e lat1, \e lon1) and the azimuth \e * azi12 are specified in the call to Rhumb::Line which returns a RhumbLine * object. RhumbLine.Position returns the location of point 2 (and, * optionally, the corresponding area, \e S12) a distance \e s12 along the * rhumb line. * * There is no public constructor for this class. (Use Rhumb::Line to create * an instance.) The Rhumb object used to create a RhumbLine must stay in * scope as long as the RhumbLine. * * Example of use: * \include example-RhumbLine.cpp **********************************************************************/ class GEOGRAPHICLIB_EXPORT RhumbLine { private: typedef Math::real real; friend class Rhumb; const Rhumb& _rh; bool _exact; real _lat1, _lon1, _azi12, _salp, _calp, _mu1, _psi1, _r1; RhumbLine& operator=(const RhumbLine&); // copy assignment not allowed RhumbLine(const Rhumb& rh, real lat1, real lon1, real azi12, bool exact); public: /** * This is a duplication of Rhumb::mask. **********************************************************************/ enum mask { /** * No output. * @hideinitializer **********************************************************************/ NONE = Rhumb::NONE, /** * Calculate latitude \e lat2. * @hideinitializer **********************************************************************/ LATITUDE = Rhumb::LATITUDE, /** * Calculate longitude \e lon2. * @hideinitializer **********************************************************************/ LONGITUDE = Rhumb::LONGITUDE, /** * Calculate azimuth \e azi12. * @hideinitializer **********************************************************************/ AZIMUTH = Rhumb::AZIMUTH, /** * Calculate distance \e s12. * @hideinitializer **********************************************************************/ DISTANCE = Rhumb::DISTANCE, /** * Calculate area \e S12. * @hideinitializer **********************************************************************/ AREA = Rhumb::AREA, /** * Unroll \e lon2 in the direct calculation. * @hideinitializer **********************************************************************/ LONG_UNROLL = Rhumb::LONG_UNROLL, /** * Calculate everything. (LONG_UNROLL is not included in this mask.) * @hideinitializer **********************************************************************/ ALL = Rhumb::ALL, }; /** * Compute the position of point 2 which is a distance \e s12 (meters) from * point 1. The area is also computed. * * @param[in] s12 distance between point 1 and point 2 (meters); it can be * negative. * @param[out] lat2 latitude of point 2 (degrees). * @param[out] lon2 longitude of point 2 (degrees). * @param[out] S12 area under the rhumb line (meters²). * * The value of \e lon2 returned is in the range [−180°, * 180°). * * If \e s12 is large enough that the rhumb line crosses a pole, the * longitude of point 2 is indeterminate (a NaN is returned for \e lon2 and * \e S12). **********************************************************************/ void Position(real s12, real& lat2, real& lon2, real& S12) const { GenPosition(s12, LATITUDE | LONGITUDE | AREA, lat2, lon2, S12); } /** * Compute the position of point 2 which is a distance \e s12 (meters) from * point 1. The area is not computed. **********************************************************************/ void Position(real s12, real& lat2, real& lon2) const { real t; GenPosition(s12, LATITUDE | LONGITUDE, lat2, lon2, t); } /** * The general position routine. RhumbLine::Position is defined in term so * this function. * * @param[in] s12 distance between point 1 and point 2 (meters); it can be * negative. * @param[in] outmask a bitor'ed combination of RhumbLine::mask values * specifying which of the following parameters should be set. * @param[out] lat2 latitude of point 2 (degrees). * @param[out] lon2 longitude of point 2 (degrees). * @param[out] S12 area under the rhumb line (meters²). * * The RhumbLine::mask values possible for \e outmask are * - \e outmask |= RhumbLine::LATITUDE for the latitude \e lat2; * - \e outmask |= RhumbLine::LONGITUDE for the latitude \e lon2; * - \e outmask |= RhumbLine::AREA for the area \e S12; * - \e outmask |= RhumbLine::ALL for all of the above; * - \e outmask |= RhumbLine::LONG_UNROLL to unroll \e lon2 instead of * wrapping it into the range [−180°, 180°). * . * With the RhumbLine::LONG_UNROLL bit set, the quantity \e lon2 − \e * lon1 indicates how many times and in what sense the rhumb line encircles * the ellipsoid. * * If \e s12 is large enough that the rhumb line crosses a pole, the * longitude of point 2 is indeterminate (a NaN is returned for \e lon2 and * \e S12). **********************************************************************/ void GenPosition(real s12, unsigned outmask, real& lat2, real& lon2, real& S12) const; /** \name Inspector functions **********************************************************************/ ///@{ /** * @return \e lat1 the latitude of point 1 (degrees). **********************************************************************/ Math::real Latitude() const { return _lat1; } /** * @return \e lon1 the longitude of point 1 (degrees). **********************************************************************/ Math::real Longitude() const { return _lon1; } /** * @return \e azi12 the azimuth of the rhumb line (degrees). **********************************************************************/ Math::real Azimuth() const { return _azi12; } /** * @return \e a the equatorial radius of the ellipsoid (meters). This is * the value inherited from the Rhumb object used in the constructor. **********************************************************************/ Math::real MajorRadius() const { return _rh.MajorRadius(); } /** * @return \e f the flattening of the ellipsoid. This is the value * inherited from the Rhumb object used in the constructor. **********************************************************************/ Math::real Flattening() const { return _rh.Flattening(); } }; } // namespace GeographicLib #endif // GEOGRAPHICLIB_RHUMB_HPP