#ifndef INTERPOLATE_POLYNOMIAL_HPP_
#define INTERPOLATE_POLYNOMIAL_HPP_

#include <vector>

#include "base/index.hpp"
#include "base/vector.hpp"

namespace VMG
{

class Grid;

class InterpolatePolynomial
{
public:
  InterpolatePolynomial(const unsigned int& degree);
  ~InterpolatePolynomial();

  void ComputeCoefficients(const Grid& grid, const Index& index);

  vmg_float Evaluate(const Vector& pos);
  void Evaluate(const Vector& pos, vmg_float& pot, Vector& field);
  Vector EvaluateNegGradient(const Vector& pos);

private:
  vmg_float& _access_coeff(const Index& index)
  {
    return coeff[index.Z() + deg_1 * (index.Y() + deg_1 * index.X())];
  }

  vmg_float& _access_coeff(const int& i, const int& j, const int& k)
  {
    return coeff[k + deg_1 * (j + deg_1 * i)];
  }

  void _compute_coefficients_1d(const Index& index, const unsigned int& direction);

  vmg_float* coeff;
  vmg_float* coeff_buffer;
  unsigned int deg, deg_1;
  Vector pos_begin;
  Vector h;
  std::vector<Vector> buffer;
  std::vector< std::vector<Vector> > buffer_diff;
};

}

#endif /* INTERPOLATE_POLYNOMIAL_HPP_ */