/**
 * @file   bspline.hpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Nov 21 13:27:22 2011
 *
 * @brief  B-Splines for molecular dynamics.
 *
 */

#ifndef BSPLINE_HPP_
#define BSPLINE_HPP_

#include "base/helper.hpp"
#include "base/index.hpp"
#include "base/particle.hpp"
#include "base/polynomial.hpp"
#include "base/vector.hpp"
#include "grid/grid.hpp"

namespace VMG
{

class BSpline
{
public:
  BSpline(const vmg_float& width);

  vmg_float EvaluateSpline(const vmg_float& val)
  {
    for (unsigned int i=0; i<intervals.size(); ++i)
      if (val < intervals[i])
	return spline_nom[i](val) / spline_denom[i](val);
    return 0.0;
  }

  void SetSpline(Grid& grid, const Particle::Particle p, const int& near_field_cells)
  {
    Index i, index_global, index_local;
    vmg_float temp_val, length;
    std::vector<vmg_float> vals;
    std::vector<vmg_float>::const_iterator iter;

    vmg_float int_val = 0.0;

    index_global = static_cast<Index>((p.Pos() - grid.Extent().Begin()) / grid.Extent().MeshWidth());
    index_local = index_global - grid.Global().BeginLocal() + grid.Local().Begin();

    // Iterate over all grid points which lie in the support of the interpolating B-Spline
    for (i.X()=-1*near_field_cells-1; i.X()<=near_field_cells+1; ++i.X())
      for (i.Y()=-1*near_field_cells-1; i.Y()<=near_field_cells+1; ++i.Y())
	for (i.Z()=-1*near_field_cells-1; i.Z()<=near_field_cells+1; ++i.Z()) {

	  // Compute distance from grid point to particle
	  length = ( p.Pos() - grid.Extent().Begin() - grid.Extent().MeshWidth() * (index_global+i) ).Length();

	  temp_val = EvaluateSpline(length);
	  vals.push_back(temp_val * p.Charge());
	  int_val += temp_val;

	}

    int_val = 1.0 / (Helper::pow(grid.MeshWidth(), 3) * int_val);

    iter = vals.begin();

    for (i.X()=-1*near_field_cells-1; i.X()<=near_field_cells+1; ++i.X())
      for (i.Y()=-1*near_field_cells-1; i.Y()<=near_field_cells+1; ++i.Y())
	for (i.Z()=-1*near_field_cells-1; i.Z()<=near_field_cells+1; ++i.Z())
	  grid(index_local+i) += *iter++ * int_val;

  }

  vmg_float EvaluatePotential(const vmg_float& val)
  {
    for (unsigned int i=0; i<intervals.size(); ++i)
      if (val < intervals[i])
	return potential_nom[i](val) / potential_denom[i](val);
    return potential_nom.back()(val) / potential_denom.back()(val);
  }

  vmg_float EvaluateDerivative(const Vector& v, const int& dv);

  const vmg_float& GetAntiDerivativeZero() const
  {
    return antid;
  }

private:
  std::vector<Polynomial> spline_nom, spline_denom;
  std::vector<Polynomial> potential_nom, potential_denom;
  vmg_float antid;
  std::vector<vmg_float> intervals;

  vmg_float R;
};

}

#endif /* BSPLINE_HPP_ */