AlbersEqualArea.hpp

Go to the documentation of this file.

 
#if !defined(GEOGRAPHICLIB_ALBERSEQUALAREA_HPP)
#define GEOGRAPHICLIB_ALBERSEQUALAREA_HPP 1
 
#include <GeographicLib/Constants.hpp>
 
namespace GeographicLib {
 
  class GEOGRAPHICLIB_EXPORT AlbersEqualArea {
  private:
    typedef Math::real real;
    real eps_, epsx_, epsx2_, tol_, tol0_;
    real _a, _f, _fm, _e2, _e, _e2m, _qZ, _qx;
    real _sign, _lat0, _k0;
    real _n0, _m02, _nrho0, _k2, _txi0, _scxi0, _sxi0;
    static const int numit_ = 5;   // Newton iterations in Reverse
    static const int numit0_ = 20; // Newton iterations in Init
    static real hyp(real x) { return Math::hypot(real(1), x); }
    // atanh(      e   * x)/      e   if f > 0
    // atan (sqrt(-e2) * x)/sqrt(-e2) if f < 0
    // x                              if f = 0
    real atanhee(real x) const {
      using std::atan2; using std::abs;
      return _f > 0 ? Math::atanh(_e * x)/_e :
        // We only invoke atanhee in txif for positive latitude.  Then x is
        // only negative for very prolate ellipsoids (_b/_a >= sqrt(2)) and we
        // still need to return a positive result in this case; hence the need
        // for the call to atan2.
        (_f < 0 ? (atan2(_e * abs(x), real(x < 0 ? -1 : 1))/_e) : x);
    }
    // return atanh(sqrt(x))/sqrt(x) - 1, accurate for small x
    static real atanhxm1(real 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
    //
    // 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);
    }
    // Datanhee(x,y) = atanhee((x-y)/(1-e^2*x*y))/(x-y)
    real Datanhee(real x, real y) const {
      real t = x - y, d = 1 - _e2 * x * y;
      return t != 0 ? atanhee(t / d) / t : 1 / d;
    }
    // DDatanhee(x,y) = (Datanhee(1,y) - Datanhee(1,x))/(y-x)
    real DDatanhee(real x, real y) const;
    void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1);
    real txif(real tphi) const;
    real tphif(real txi) const;
 
    friend class Ellipsoid;           // For access to txif, tphif, etc.
  public:
 
    AlbersEqualArea(real a, real f, real stdlat, real k0);
 
    AlbersEqualArea(real a, real f, real stdlat1, real stdlat2, real k1);
 
    AlbersEqualArea(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 AlbersEqualArea& CylindricalEqualArea();
 
    static const AlbersEqualArea& AzimuthalEqualAreaNorth();
 
    static const AlbersEqualArea& AzimuthalEqualAreaSouth();
  };
 
} // namespace GeographicLib
 
#endif  // GEOGRAPHICLIB_ALBERSEQUALAREA_HPP