Rhumb.hpp

Go to the documentation of this file.

#if !defined(GEOGRAPHICLIB_RHUMB_HPP)
#define GEOGRAPHICLIB_RHUMB_HPP 1

#include <GeographicLib/Constants.hpp>
#include <GeographicLib/Ellipsoid.hpp>

#if !defined(GEOGRAPHICLIB_RHUMBAREA_ORDER)

#  define GEOGRAPHICLIB_RHUMBAREA_ORDER \
  (GEOGRAPHICLIB_PRECISION == 2 ? 6 : \
   (GEOGRAPHICLIB_PRECISION == 1 ? 4 : 8))
#endif

namespace GeographicLib {

  class RhumbLine;
  template <class T> class PolygonAreaT;

  class  GEOGRAPHICLIB_EXPORT Rhumb {
  private:
    typedef Math::real real;
    friend class RhumbLine;
    template <class T> 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 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://doi.org/10.1145/334714.334716
    //   http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf

    static real Dlog(real x, real y) {
      real t = x - y;
      return t != 0 ? 2 * Math::atanh(t / (x + y)) / t : 1 / x;
    }
    // N.B., x and y are in degrees
    static real Dtan(real x, real y) {
      real d = x - y, tx = Math::tand(x), ty = Math::tand(y), txy = tx * ty;
      return d != 0 ?
        (2 * txy > -1 ? (1 + txy) * Math::tand(d) : tx - ty) /
        (d * Math::degree()) :
        1 + txy;
    }
    static real Datan(real x, real y) {
      using std::atan;
      real d = x - y, xy = x * y;
      return d != 0 ?
        (2 * xy > -1 ? atan( d / (1 + xy) ) : atan(x) - atan(y)) / d :
        1 / (1 + xy);
    }
    static real Dsin(real x, real y) {
      using std::sin; using std::cos;
      real d = (x - y) / 2;
      return cos((x + y)/2) * (d != 0 ? sin(d) / d : 1);
    }
    static real Dsinh(real x, real y) {
      using std::sinh; using std::cosh;
      real d = (x - y) / 2;
      return cosh((x + y) / 2) * (d != 0 ? sinh(d) / d : 1);
    }
    static real Dcosh(real x, real y) {
      using std::sinh;
      real d = (x - y) / 2;
      return sinh((x + y) / 2) * (d != 0 ? sinh(d) / d : 1);
    }
    static real Dasinh(real x, real y) {
      real d = x - y,
        hx = Math::hypot(real(1), x), hy = Math::hypot(real(1), y);
      return d != 0 ? Math::asinh(x*y > 0 ? d * (x + y) / (x*hy + y*hx) :
                                  x*hy - y*hx) / d :
        1 / hx;
    }
    static 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 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)
    real Deatanhe(real x, real y) const {
      real t = x - y, d = 1 - _ell._e2 * x * y;
      return t != 0 ? 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:

    enum mask {
      NONE          = 0U,
      LATITUDE      = 1U<<7,
      LONGITUDE     = 1U<<8,
      AZIMUTH       = 1U<<9,
      DISTANCE      = 1U<<10,
      AREA          = 1U<<14,
      LONG_UNROLL   = 1U<<15,
      ALL           = 0x7F80U,
    };

    Rhumb(real a, real f, bool exact = true);

    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);
    }

    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);
    }

    void GenDirect(real lat1, real lon1, real azi12, real s12,
                   unsigned outmask, real& lat2, real& lon2, real& S12) const;

    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);
    }

    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);
    }

    void GenInverse(real lat1, real lon1, real lat2, real lon2,
                    unsigned outmask,
                    real& s12, real& azi12, real& S12) const;

    RhumbLine Line(real lat1, real lon1, real azi12) const;


    Math::real MajorRadius() const { return _ell.MajorRadius(); }

    Math::real Flattening() const { return _ell.Flattening(); }

    Math::real EllipsoidArea() const { return _ell.Area(); }

    static const Rhumb& WGS84();
  };

  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:

    enum mask {
      NONE          = Rhumb::NONE,
      LATITUDE      = Rhumb::LATITUDE,
      LONGITUDE     = Rhumb::LONGITUDE,
      AZIMUTH       = Rhumb::AZIMUTH,
      DISTANCE      = Rhumb::DISTANCE,
      AREA          = Rhumb::AREA,
      LONG_UNROLL   = Rhumb::LONG_UNROLL,
      ALL           = Rhumb::ALL,
    };

    void Position(real s12, real& lat2, real& lon2, real& S12) const {
      GenPosition(s12, LATITUDE | LONGITUDE | AREA, lat2, lon2, S12);
    }

    void Position(real s12, real& lat2, real& lon2) const {
      real t;
      GenPosition(s12, LATITUDE | LONGITUDE, lat2, lon2, t);
    }

    void GenPosition(real s12, unsigned outmask,
                     real& lat2, real& lon2, real& S12) const;


    Math::real Latitude() const { return _lat1; }

    Math::real Longitude() const { return _lon1; }

    Math::real Azimuth() const { return  _azi12; }

    Math::real MajorRadius() const { return _rh.MajorRadius(); }

    Math::real Flattening() const { return _rh.Flattening(); }
  };

} // namespace GeographicLib

#endif  // GEOGRAPHICLIB_RHUMB_HPP