////////////////////////////////////////////////////////////////////
/*
MultivariateInterp(npoints, vector x, vector y, vector z) approximates the value f(x,y) = z via inverse distance weighting based on the LENGTH-1 nearest data points on an irregular grid (modified Shepard's method using only nearest data points within sphere with radius=LENGTH).
Could be made faster by using kdtree.
*/

#ifndef __RAT_MultivariateInterp__
#define __RAT_MultivariateInterp__

#include <vector>
#define LENGTH 5

namespace RAT {


template <class NumType>
class MultivariateInterp {
public:
  MultivariateInterp() { npoints = 0; nneighbours = 0; };
  MultivariateInterp(int _npoints, int _nneighbours, const NumType _x[], const NumType _y[], const NumType _z[]) {
    Set(_npoints, _nneighbours, _x, _y, _z);
  };
  MultivariateInterp(const std::vector<NumType> &_x, const std::vector<NumType> &_y, const std::vector<NumType> &_z) {
    Set(_x, _y, _z);
  };

  void Set(int _npoints, int _nneighbours, const NumType _x[], const NumType _y[], const NumType _z[]);
  void Set(const std::vector<NumType> &_x, const std::vector<NumType> &_y, const std::vector<NumType> &_z) {
    npoints = _x.size(); nneighbours = LENGTH; x = _x; y = _y; z = _z;
  };

  NumType operator() (const NumType xeval, const NumType yeval);

  int Points() { return npoints; };

protected:
  int* GetNN(const NumType xeval, const NumType yeval);

  int npoints, nneighbours;
  std::vector<NumType> x, y, z;
};

#include "MultivariateInterp.icc"

} // namespace RAT

#endif