LambertConformalConic.hpp

Go to the documentation of this file.

#if !defined(GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP)
#define GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP 1

#include <GeographicLib/Constants.hpp>

namespace GeographicLib {

  class GEOGRAPHICLIB_EXPORT LambertConformalConic {
  private:
    typedef Math::real real;
    real eps_, epsx_, ahypover_;
    real _a, _f, _fm, _e2, _es;
    real _sign, _n, _nc, _t0nm1, _scale, _lat0, _k0;
    real _scbet0, _tchi0, _scchi0, _psi0, _nrho0, _drhomax;
    static const int numit_ = 5;
    static real hyp(real x) { return Math::hypot(real(1), x); }
    // Divided differences
    // Definition: Df(x,y) = (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
    //
    // General rules
    // h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y)
    // h(x) = f(x)*g(x):
    //        Dh(x,y) = Df(x,y)*g(x) + Dg(x,y)*f(y)
    //                = Df(x,y)*g(y) + Dg(x,y)*f(x)
    //                = Df(x,y)*(g(x)+g(y))/2 + Dg(x,y)*(f(x)+f(y))/2
    //
    // hyp(x) = sqrt(1+x^2): Dhyp(x,y) = (x+y)/(hyp(x)+hyp(y))
    static real Dhyp(real x, real y, real hx, real hy)
    // hx = hyp(x)
    { return (x + y) / (hx + hy); }
    // sn(x) = x/sqrt(1+x^2): Dsn(x,y) = (x+y)/((sn(x)+sn(y))*(1+x^2)*(1+y^2))
    static real Dsn(real x, real y, real sx, real sy) {
      // sx = x/hyp(x)
      real t = x * y;
      return t > 0 ? (x + y) * Math::sq( (sx * sy)/t ) / (sx + sy) :
        (x - y != 0 ? (sx - sy) / (x - y) : 1);
    }
    // Dlog1p(x,y) = log1p((x-y)/(1+y))/(x-y)
    static real Dlog1p(real x, real y) {
      real t = x - y; if (t < 0) { t = -t; y = x; }
      return t != 0 ? Math::log1p(t / (1 + y)) / t : 1 / (1 + x);
    }
    // Dexp(x,y) = exp((x+y)/2) * 2*sinh((x-y)/2)/(x-y)
    static real Dexp(real x, real y) {
      using std::sinh; using std::exp;
      real t = (x - y)/2;
      return (t != 0 ? sinh(t)/t : 1) * exp((x + y)/2);
    }
    // Dsinh(x,y) = 2*sinh((x-y)/2)/(x-y) * cosh((x+y)/2)
    //   cosh((x+y)/2) = (c+sinh(x)*sinh(y)/c)/2
    //   c=sqrt((1+cosh(x))*(1+cosh(y)))
    //   cosh((x+y)/2) = sqrt( (sinh(x)*sinh(y) + cosh(x)*cosh(y) + 1)/2 )
    static real Dsinh(real x, real y, real sx, real sy, real cx, real cy)
    // sx = sinh(x), cx = cosh(x)
    {
      // real t = (x - y)/2, c = sqrt((1 + cx) * (1 + cy));
      // return (t ? sinh(t)/t : real(1)) * (c + sx * sy / c) /2;
      using std::sinh; using std::sqrt;
      real t = (x - y)/2;
      return (t != 0 ? sinh(t)/t : 1) * sqrt((sx * sy + cx * cy + 1) /2);
    }
    // Dasinh(x,y) = asinh((x-y)*(x+y)/(x*sqrt(1+y^2)+y*sqrt(1+x^2)))/(x-y)
    //             = asinh((x*sqrt(1+y^2)-y*sqrt(1+x^2)))/(x-y)
    static real Dasinh(real x, real y, real hx, real hy) {
      // hx = hyp(x)
      real t = x - y;
      return t != 0 ?
        Math::asinh(x*y > 0 ? t * (x+y) / (x*hy + y*hx) : x*hy - y*hx) / t :
        1/hx;
    }
    // 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 - _e2 * x * y;
      return t != 0 ? Math::eatanhe(t / d, _es) / t : _e2 / d;
    }
    void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1);
  public:

    LambertConformalConic(real a, real f, real stdlat, real k0);

    LambertConformalConic(real a, real f, real stdlat1, real stdlat2, real k1);

    LambertConformalConic(real a, real f,
                          real sinlat1, real coslat1,
                          real sinlat2, real coslat2,
                          real k1);

    void SetScale(real lat, real k = real(1));

    void Forward(real lon0, real lat, real lon,
                 real& x, real& y, real& gamma, real& k) const;

    void Reverse(real lon0, real x, real y,
                 real& lat, real& lon, real& gamma, real& k) const;

    void Forward(real lon0, real lat, real lon,
                 real& x, real& y) const {
      real gamma, k;
      Forward(lon0, lat, lon, x, y, gamma, k);
    }

    void Reverse(real lon0, real x, real y,
                 real& lat, real& lon) const {
      real gamma, k;
      Reverse(lon0, x, y, lat, lon, gamma, k);
    }


    Math::real MajorRadius() const { return _a; }

    Math::real Flattening() const { return _f; }

    Math::real OriginLatitude() const { return _lat0; }

    Math::real CentralScale() const { return _k0; }

    static const LambertConformalConic& Mercator();
  };

} // namespace GeographicLib

#endif  // GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP