SphericalEngine.hpp

Go to the documentation of this file.

#if !defined(GEOGRAPHICLIB_SPHERICALENGINE_HPP)
#define GEOGRAPHICLIB_SPHERICALENGINE_HPP 1

#include <vector>
#include <istream>
#include <GeographicLib/Constants.hpp>

#if defined(_MSC_VER)
// Squelch warnings about dll vs vector
#  pragma warning (push)
#  pragma warning (disable: 4251)
#endif

namespace GeographicLib {

  class CircularEngine;

  class GEOGRAPHICLIB_EXPORT SphericalEngine {
  private:
    typedef Math::real real;
    // CircularEngine needs access to sqrttable, scale
    friend class CircularEngine;
    // Return the table of the square roots of integers
    static std::vector<real>& sqrttable();
    // An internal scaling of the coefficients to avoid overflow in
    // intermediate calculations.
    static real scale() {
      using std::pow;
      static const real
        // Need extra real because, since C++11, pow(float, int) returns double
        s = real(pow(real(std::numeric_limits<real>::radix),
                     -3 * (std::numeric_limits<real>::max_exponent < (1<<14) ?
                           std::numeric_limits<real>::max_exponent : (1<<14))
                     / 5));
      return s;
    }
    // Move latitudes near the pole off the axis by this amount.
    static real eps() {
      using std::sqrt;
      return std::numeric_limits<real>::epsilon() *
        sqrt(std::numeric_limits<real>::epsilon());
    }
    SphericalEngine();          // Disable constructor
  public:
    enum normalization {
      FULL = 0,
      SCHMIDT = 1,
    };

    class GEOGRAPHICLIB_EXPORT coeff {
    private:
      int _Nx, _nmx, _mmx;
      std::vector<real>::const_iterator _Cnm;
      std::vector<real>::const_iterator _Snm;
    public:
      coeff() : _Nx(-1) , _nmx(-1) , _mmx(-1) {}
      coeff(const std::vector<real>& C,
            const std::vector<real>& S,
            int N, int nmx, int mmx)
        : _Nx(N)
        , _nmx(nmx)
        , _mmx(mmx)
        , _Cnm(C.begin())
        , _Snm(S.begin())
      {
        if (!(_Nx >= _nmx && _nmx >= _mmx && _mmx >= -1))
          throw GeographicErr("Bad indices for coeff");
        if (!(index(_nmx, _mmx) < int(C.size()) &&
              index(_nmx, _mmx) < int(S.size()) + (_Nx + 1)))
          throw GeographicErr("Arrays too small in coeff");
        SphericalEngine::RootTable(_nmx);
      }
      coeff(const std::vector<real>& C,
            const std::vector<real>& S,
            int N)
        : _Nx(N)
        , _nmx(N)
        , _mmx(N)
        , _Cnm(C.begin())
        , _Snm(S.begin())
      {
        if (!(_Nx >= -1))
          throw GeographicErr("Bad indices for coeff");
        if (!(index(_nmx, _mmx) < int(C.size()) &&
              index(_nmx, _mmx) < int(S.size()) + (_Nx + 1)))
          throw GeographicErr("Arrays too small in coeff");
        SphericalEngine::RootTable(_nmx);
      }
      int N() const { return _Nx; }
      int nmx() const { return _nmx; }
      int mmx() const { return _mmx; }
      int index(int n, int m) const
      { return m * _Nx - m * (m - 1) / 2 + n; }
      Math::real Cv(int k) const { return *(_Cnm + k); }
      Math::real Sv(int k) const { return *(_Snm + (k - (_Nx + 1))); }
      Math::real Cv(int k, int n, int m, real f) const
      { return m > _mmx || n > _nmx ? 0 : *(_Cnm + k) * f; }
      Math::real Sv(int k, int n, int m, real f) const
      { return m > _mmx || n > _nmx ? 0 : *(_Snm + (k - (_Nx + 1))) * f; }

      static int Csize(int N, int M)
      { return (M + 1) * (2 * N - M + 2) / 2; }

      static int Ssize(int N, int M)
      { return Csize(N, M) - (N + 1); }

      static void readcoeffs(std::istream& stream, int& N, int& M,
                             std::vector<real>& C, std::vector<real>& S);
    };

    template<bool gradp, normalization norm, int L>
      static Math::real Value(const coeff c[], const real f[],
                              real x, real y, real z, real a,
                              real& gradx, real& grady, real& gradz);

    template<bool gradp, normalization norm, int L>
      static CircularEngine Circle(const coeff c[], const real f[],
                                   real p, real z, real a);
    static void RootTable(int N);

    static void ClearRootTable() {
      std::vector<real> temp(0);
      sqrttable().swap(temp);
    }
  };

} // namespace GeographicLib

#if defined(_MSC_VER)
#  pragma warning (pop)
#endif

#endif  // GEOGRAPHICLIB_SPHERICALENGINE_HPP