JacobiConformal.hpp

Go to the documentation of this file.

#include <GeographicLib/EllipticFunction.hpp>

namespace GeographicLib {
  class JacobiConformal {
    typedef Math::real real;
    real _a, _b, _c, _ab2, _bc2, _ac2;
    EllipticFunction _ex, _ey;
    static void norm(real& x, real& y)
    { real z = Math::hypot(x, y); x /= z; y /= z; }
  public:
    JacobiConformal(real a, real b, real c)
      : _a(a), _b(b), _c(c)
      , _ab2((_a - _b) * (_a + _b))
      , _bc2((_b - _c) * (_b + _c))
      , _ac2((_a - _c) * (_a + _c))
      , _ex(_ab2 / _ac2 * Math::sq(_c / _b), -_ab2 / Math::sq(_b),
            _bc2 / _ac2 * Math::sq(_a / _b), Math::sq(_a / _b))
      , _ey(_bc2 / _ac2 * Math::sq(_a / _b), +_bc2 / Math::sq(_b),
            _ab2 / _ac2 * Math::sq(_c / _b), Math::sq(_c / _b))
    {
      if (!(Math::isfinite(_a) && _a >= _b && _b >= _c && _c > 0))
        throw GeographicErr("JacobiConformal: axes are not in order");
      if (!(_a > _c))
        throw GeographicErr
          ("JacobiConformal: use alternate constructor for sphere");
    }
    JacobiConformal(real a, real b, real c, real ab, real bc)
      : _a(a), _b(b), _c(c)
      , _ab2(ab * (_a + _b))
      , _bc2(bc * (_b + _c))
      , _ac2(_ab2 + _bc2)
      , _ex(_ab2 / _ac2 * Math::sq(_c / _b),
            -(_a - _b) * (_a + _b) / Math::sq(_b),
            _bc2 / _ac2 * Math::sq(_a / _b), Math::sq(_a / _b))
      , _ey(_bc2 / _ac2 * Math::sq(_a / _b),
            +(_b - _c) * (_b + _c) / Math::sq(_b),
            _ab2 / _ac2 * Math::sq(_c / _b), Math::sq(_c / _b))
    {
      if (!(Math::isfinite(_a) && _a >= _b && _b >= _c && _c > 0 &&
            ab >= 0 && bc >= 0))
        throw GeographicErr("JacobiConformal: axes are not in order");
      if (!(ab + bc > 0 && Math::isfinite(_ac2)))
        throw GeographicErr("JacobiConformal: ab + bc must be positive");
    }
    Math::real x() const { return Math::sq(_a / _b) * _ex.Pi(); }
    Math::real x(real somg, real comg) const {
      real somg1 = _b * somg, comg1 = _a * comg; norm(somg1, comg1);
      return Math::sq(_a / _b)
        * _ex.Pi(somg1, comg1, _ex.Delta(somg1, comg1));
    }
    Math::real x(real omg) const {
      real somg, comg;
      Math::sincosd(omg, somg, comg);
      return x(somg, comg) / Math::degree();
    }
    Math::real y() const { return Math::sq(_c / _b) * _ey.Pi(); }
    Math::real y(real sbet, real cbet) const {
      real sbet1 = _b * sbet, cbet1 = _c * cbet; norm(sbet1, cbet1);
      return Math::sq(_c / _b)
        * _ey.Pi(sbet1, cbet1, _ey.Delta(sbet1, cbet1));
    }
    Math::real y(real bet) const {
      real sbet, cbet;
      Math::sincosd(bet, sbet, cbet);
      return y(sbet, cbet) / Math::degree();
    }
  };

} // namespace GeographicLib