diff --git a/src/interpolation/interpolation.hpp b/src/interpolation/interpolation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7af561b114dd793a49c4fa09f81cabd456de8766
--- /dev/null
+++ b/src/interpolation/interpolation.hpp
@@ -0,0 +1,367 @@
+/*
+ * mp4_interpolation.hpp
+ *
+ * Created on: May 4, 2017
+ * Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_HPP_
+#define OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_HPP_
+
+#include "NN/CellList/MemFast.hpp"
+#include "NN/CellList/CellList.hpp"
+
+#define INTERPOLATION_ERROR_OBJECT std::runtime_error("Runtime interpolation error");
+
+constexpr int inte_m2p = 0;
+constexpr int inte_p2m = 1;
+
+template<typename T>
+struct mul_inte
+{
+ inline static void value(T & result, const T & coeff, const T & src)
+ {
+ result += coeff * src;
+ }
+};
+
+template<typename T, unsigned int N1>
+struct mul_inte<T[N1]>
+{
+ inline static void value(T (& result)[N1], const T & coeff, const T (& src)[N1])
+ {
+ for (size_t i = 0 ; i < N1 ; i++)
+ result[i] += coeff * src[i];
+ }
+};
+
+template<typename T, unsigned int N1, unsigned int N2>
+struct mul_inte<T[N1][N2]>
+{
+ inline static void value(T (& result)[N1][N2], const T & coeff, const T (& src)[N1][N2])
+ {
+ for (size_t i = 0 ; i < N1 ; i++)
+ for (size_t j = 0 ; j < N2 ; j++)
+ result[i][j] += coeff * src[i][j];
+ }
+};
+
+template<unsigned int prp_g, unsigned int prp_v,unsigned int m2p_or_p2m>
+struct inte_template
+{
+ template<typename grid, typename vector, typename iterator, typename key_type> inline static void value(grid & gd,
+ vector & vd,
+ const grid_dist_key_dx<vector::dims> & k_dist,
+ iterator & key_p,
+ const grid_cpu<vector::dims,aggregate<typename vector::stype>> & a_int,
+ const key_type & key)
+ {
+ mul_inte<typename std::remove_reference<decltype(gd.template get<prp_g>(k_dist))>::type>::value(gd.template get<prp_g>(k_dist),a_int.template get<0>(key),vd.template getProp<prp_v>(key_p));
+ }
+};
+
+template<unsigned int prp_g, unsigned int prp_v>
+struct inte_template<prp_g,prp_v,inte_m2p>
+{
+ template<typename grid, typename vector, typename iterator, typename key_type> inline static void value(grid & gd,
+ vector & vd,
+ const grid_dist_key_dx<vector::dims> & k_dist,
+ iterator & key_p,
+ const grid_cpu<vector::dims,aggregate<typename vector::stype>> & a_int,
+ const key_type & key)
+ {
+ mul_inte<typename std::remove_reference<decltype(gd.template get<prp_g>(k_dist))>::type>::value(vd.template getProp<prp_v>(key_p),a_int.template get<0>(key),gd.template get<prp_g>(k_dist));
+ }
+};
+
+
+/*! \brief return the sub-domain where this particle must be interpolated
+ *
+ * \param p particle position
+ *
+ * \return the sub-domain id
+ *
+ */
+template<typename vector, typename grid>
+inline size_t getSub(Point<vector::dims,typename vector::stype> & p,
+ const CellList<vector::dims,typename vector::stype,Mem_fast<vector::dims,typename vector::stype>,shift<vector::dims,typename vector::stype>> & geo_cell,
+ grid & gd)
+{
+ size_t cell = geo_cell.getCell(p);
+
+ for (size_t i = 0 ; i < geo_cell.getNelements(cell) ; i++)
+ {
+ size_t ns = geo_cell.get(cell,i);
+
+ if (gd.getDecomposition().getSubDomain(ns).isInside(p))
+ return ns;
+ }
+
+#ifdef SE_CLASS1
+ std::cerr << __FILE__ << ":" << __LINE__ << " Error: " << std::endl;
+ ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
+#endif
+
+ return (size_t)-1;
+}
+
+
+template<unsigned int prp_g, unsigned int prp_v, unsigned int m2p_or_p2m, typename kernel, typename iterator, typename vector, typename grid, typename grid_inte>
+ inline void inte_calc(iterator & it,
+ const vector & vd,
+ const Box<vector::dims,typename vector::stype> & domain,
+ int (& ip)[vector::dims][kernel::np],
+ grid & gd,
+ const typename vector::stype (& dx)[vector::dims],
+ typename vector::stype (& xp)[vector::dims],
+ const grid_inte & a_int,
+ typename vector::stype (& a)[vector::dims][kernel::np],
+ typename vector::stype (& x)[vector::dims][kernel::np],
+ size_t (& sz)[vector::dims],
+ const CellList<vector::dims,typename vector::stype,Mem_fast<vector::dims,typename vector::stype>,shift<vector::dims,typename vector::stype>> & geo_cell)
+{
+ auto key_p = it.get();
+
+ Point<vector::dims,typename vector::stype> p = vd.getPos(key_p);
+
+ // On which sub-domain we interpolate the particle
+ size_t sub = getSub<vector>(p,geo_cell,gd);
+
+ typename vector::stype x0[vector::dims];
+
+ // calculate the position of the particle in the grid
+ // coordinated
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ x0[i] = (p.get(i)-domain.getLow(i))*dx[i];
+
+ // convert into integer
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ ip[i][0] = (int)x0[i];
+
+ // convert the global grid position into local grid position
+ grid_key_dx<vector::dims> base;
+
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ base.set_d(i,ip[i][0] - gd.getLocalGridsInfo().get(sub).origin.get(i) - (long int)kernel::np/2 + 1);
+
+ // convenient grid of points
+
+ for (size_t j = 0 ; j < kernel::np-1 ; j++)
+ {
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ ip[i][j+1] = (int)ip[i][j]+1;
+ }
+
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ xp[i] = x0[i] - ip[i][0];
+
+ for (long int j = 0 ; j < kernel::np ; j++)
+ {
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ x[i][j] = - xp[i] + typename vector::stype((long int)j - (long int)kernel::np/2 + 1);
+ }
+
+ for (size_t j = 0 ; j < kernel::np ; j++)
+ {
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ a[i][j] = kernel::value(x[i][j],j);
+ }
+
+ grid_sm<vector::dims,void> gs(sz);
+ grid_key_dx_iterator<vector::dims> kit(gs);
+
+ while (kit.isNext())
+ {
+ auto key = kit.get();
+
+ a_int.template get<0>(key) = 1;
+
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ a_int.template get<0>(key) *= a[i][key.get(i)];
+
+ ++kit;
+ }
+
+ grid_key_dx_iterator<vector::dims> kit2(gs);
+ grid_dist_key_dx<vector::dims> k_dist;
+ k_dist.setSub(sub);
+
+ while (kit2.isNext())
+ {
+ auto key = kit2.get();
+
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ k_dist.getKeyRef().set_d(i,key.get(i) + base.get(i));
+
+ inte_template<prp_g,prp_v,m2p_or_p2m>::value(gd,vd,k_dist,key_p,a_int,key);
+
+ ++kit2;
+ }
+}
+
+
+template<typename vector,typename grid, typename kernel>
+class interpolate
+{
+ CellList<vector::dims,typename vector::stype,Mem_fast<vector::dims,typename vector::stype>,shift<vector::dims,typename vector::stype>> geo_cell;
+
+ struct Box_vol
+ {
+ Box<vector::dims,size_t> bv;
+
+ size_t vol;
+
+ void operator<(const Box_vol & bv)
+ {
+ return vol < bv.vol;
+ }
+ };
+
+ //! particles
+ vector & vd;
+
+ //! grid
+ grid & gd;
+
+
+public:
+
+ interpolate(vector & vd, grid & gd)
+ :vd(vd),gd(gd)
+ {
+ // get the processor bounding box in grid units
+ Box<vector::dims,typename vector::stype> bb = gd.getDecomposition().getProcessorBounds();
+ Box<vector::dims,typename vector::stype> bunit = gd.getDecomposition().getCellDecomposer().getCellBox();
+
+ size_t div[vector::dims];
+
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ div[i] = (bb.getHigh(i) - bb.getLow(i)) / bunit.getHigh(i);
+
+ geo_cell.Initialize(bb,div);
+
+ // Now draw the domain into the cell list
+
+ auto & dec = gd.getDecomposition();
+
+ for (size_t i = 0 ; i < dec.getNSubDomain() ; i++)
+ {
+ const Box<vector::dims,typename vector::stype> & bx = dec.getSubDomain(i);
+
+ // get the cells this box span
+ const grid_key_dx<vector::dims> p1 = geo_cell.getCellGrid(bx.getP1());
+ const grid_key_dx<vector::dims> p2 = geo_cell.getCellGrid(bx.getP2());
+
+ // Get the grid and the sub-iterator
+ auto & gi = geo_cell.getGrid();
+ grid_key_dx_iterator_sub<vector::dims> g_sub(gi,p1,p2);
+
+ // add the box-id to the cell list
+ while (g_sub.isNext())
+ {
+ auto key = g_sub.get();
+ geo_cell.addCell(gi.LinId(key),i);
+ ++g_sub;
+ }
+ }
+ };
+
+ template<unsigned int prp_v, unsigned int prp_g> void p2m(vector & vd, grid & gd)
+ {
+#ifdef SE_CLASS1
+
+ if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
+ {
+ std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
+ " and the grid is different. In order to interpolate the two data structure must have the" <<
+ " same decomposition" << std::endl;
+
+ ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
+ }
+
+#endif
+
+ Box<vector::dims,typename vector::stype> domain = vd.getDecomposition().getDomain();
+
+ // grid spacing
+ typename vector::stype dx[vector::dims];
+
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ dx[i] = 1.0/gd.spacing(i);
+
+ // point position
+ typename vector::stype xp[vector::dims];
+
+ int ip[vector::dims][kernel::np];
+ typename vector::stype x[vector::dims][kernel::np];
+ typename vector::stype a[vector::dims][kernel::np];
+
+ size_t sz[vector::dims];
+
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ sz[i] = kernel::np;
+
+ grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
+ a_int.setMemory();
+
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext() == true)
+ {
+ inte_calc<prp_g,prp_v,inte_p2m,kernel>(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell);
+
+ ++it;
+ }
+ }
+
+ template<unsigned int prp_v, unsigned int prp_g> void m2p(grid & gd, vector & vd)
+ {
+#ifdef SE_CLASS1
+
+ if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
+ {
+ std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
+ " and the grid is different. In order to interpolate the two data structure must have the" <<
+ " same decomposition" << std::endl;
+
+ ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
+ }
+
+#endif
+
+ Box<vector::dims,typename vector::stype> domain = vd.getDecomposition().getDomain();
+
+ // grid spacing
+ typename vector::stype dx[vector::dims];
+
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ dx[i] = 1.0/gd.spacing(i);
+
+ // point position
+ typename vector::stype xp[vector::dims];
+
+ int ip[vector::dims][kernel::np];
+ typename vector::stype x[vector::dims][kernel::np];
+ typename vector::stype a[vector::dims][kernel::np];
+
+ size_t sz[vector::dims];
+
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ sz[i] = kernel::np;
+
+ grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
+ a_int.setMemory();
+
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext() == true)
+ {
+ inte_calc<prp_g,prp_v,inte_m2p,kernel>(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell);
+
+ ++it;
+ }
+ }
+
+};
+
+#endif /* OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_HPP_ */
diff --git a/src/interpolation/interpolation_unit_tests.hpp b/src/interpolation/interpolation_unit_tests.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..888341672b4dc113dd297af2cb4f030e0b5b014a
--- /dev/null
+++ b/src/interpolation/interpolation_unit_tests.hpp
@@ -0,0 +1,423 @@
+/*
+ * interpolation_unit_tests.hpp
+ *
+ * Created on: May 5, 2017
+ * Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_UNIT_TESTS_HPP_
+#define OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_UNIT_TESTS_HPP_
+
+#include "interpolation/mp4_kernel.hpp"
+#include "interpolation/interpolation.hpp"
+#include "interpolation/z_spline.hpp"
+
+BOOST_AUTO_TEST_SUITE( interpolation_test )
+
+template<typename grid, unsigned int mom_p> void momenta_grid(grid & gd,typename grid::stype (& mom_tot)[grid::dims])
+{
+ auto it = gd.getDomainGhostIterator();
+
+ for (size_t i = 0 ; i < grid::dims ; i++)
+ mom_tot[i] = 0.0;
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+ auto key_g = gd.getGKey(key);
+
+ for (size_t i = 0 ; i < grid::dims ; i++)
+ {
+ typename grid::stype coord = gd.spacing(i)*key_g.get(i);
+
+ mom_tot[i] += boost::math::pow<mom_p>(coord)*gd.template getProp<0>(key);
+ }
+
+ ++it;
+ }
+}
+
+template<typename grid, unsigned int mom_p> void momenta_grid_domain(grid & gd,typename grid::stype (& mom_tot)[grid::dims])
+{
+ auto it = gd.getDomainIterator();
+
+ for (size_t i = 0 ; i < grid::dims ; i++)
+ mom_tot[i] = 0.0;
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+ auto key_g = gd.getGKey(key);
+
+ for (size_t i = 0 ; i < grid::dims ; i++)
+ {
+ typename grid::stype coord = gd.spacing(i)*key_g.get(i);
+
+ mom_tot[i] += boost::math::pow<mom_p>(coord)*gd.template getProp<0>(key);
+ }
+
+ ++it;
+ }
+}
+
+template<typename vector, unsigned int mom_p> void momenta_vector(vector & vd,typename vector::stype (& mom_tot)[vector::dims])
+{
+ auto it = vd.getDomainIterator();
+
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ mom_tot[i] = 0.0;
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ for (size_t i = 0 ; i < vector::dims ; i++)
+ {
+ typename vector::stype coord = vd.getPos(key)[i];
+
+ mom_tot[i] += boost::math::pow<mom_p>(coord)*vd.template getProp<0>(key);
+ }
+
+ ++it;
+ }
+}
+
+
+BOOST_AUTO_TEST_CASE( interpolation_full_test )
+{
+ Box<2,float> domain({0.0,0.0},{1.0,1.0});
+ size_t sz[2] = {64,64};
+
+ Ghost<2,long int> gg(3);
+ Ghost<2,float> gv(0.01);
+
+ size_t bc_v[2] = {PERIODIC,PERIODIC};
+
+ {
+ vector_dist<2,float,aggregate<float>> vd(4096,domain,bc_v,gv);
+ grid_dist_id<2,float,aggregate<float>> gd(vd.getDecomposition(),sz,gg);
+
+ // set one particle on vd
+
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto p = it.get();
+
+ vd.getPos(p)[0] = (double)rand()/RAND_MAX;
+ vd.getPos(p)[1] = (double)rand()/RAND_MAX;
+ vd.getPos(p)[2] = (double)rand()/RAND_MAX;
+
+ vd.getProp<0>(p) = 5.0/*(double)rand()/RAND_MAX*/;
+
+ ++it;
+ }
+
+ vd.map();
+
+ // Reset the grid
+
+ auto it2 = gd.getDomainGhostIterator();
+
+ while (it2.isNext())
+ {
+ auto key = it2.get();
+
+ gd.template get<0>(key) = 0.0;
+
+ ++it2;
+ }
+
+ interpolate<decltype(vd),decltype(gd),mp4_kernel<float>> inte(vd,gd);
+
+ inte.p2m<0,0>(vd,gd);
+
+ float mg[2];
+ float mv[2];
+
+ momenta_grid<decltype(gd),0>(gd,mg);
+ momenta_vector<decltype(vd),0>(vd,mv);
+
+ BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
+ BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+
+ momenta_grid<decltype(gd),1>(gd,mg);
+ momenta_vector<decltype(vd),1>(vd,mv);
+
+ BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
+ BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+
+ momenta_grid<decltype(gd),2>(gd,mg);
+ momenta_vector<decltype(vd),2>(vd,mv);
+
+ BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
+ BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+
+ // We have to do a ghost get before interpolating m2p
+ // Before doing mesh to particle particle must be arranged
+ // into a grid like
+
+ vd.clear();
+
+ auto it4 = vd.getGridIterator(sz);
+
+ while(it4.isNext())
+ {
+ auto key = it4.get();
+
+ vd.add();
+
+ vd.getLastPos()[0] = key.get(0) * it4.getSpacing(0) + domain.getLow(0) + 0.1*it4.getSpacing(0);
+ vd.getLastPos()[1] = key.get(1) * it4.getSpacing(1) + domain.getLow(1) + 0.1*it4.getSpacing(1);
+
+ vd.getLastProp<0>() = 0.0;
+
+ ++it4;
+ }
+
+ // Reset also the grid
+
+ auto it5 = gd.getDomainGhostIterator();
+
+ while(it5.isNext())
+ {
+ auto key = it5.get();
+
+ gd.get<0>(key) = 0.0;
+
+ ++it5;
+ }
+ gd.ghost_get<0>();
+
+ grid_key_dx<2> start({3,3});
+ grid_key_dx<2> stop({(long int)gd.size(0) - 3,(long int)gd.size(1) - 3});
+
+ auto it6 = gd.getSubDomainIterator(start,stop);
+ while(it6.isNext())
+ {
+ auto key = it6.get();
+
+ gd.get<0>(key) = (double)rand()/RAND_MAX;;
+
+ ++it6;
+ }
+ gd.ghost_get<0>();
+
+ vd.map();
+ gd.ghost_get<0>();
+ inte.m2p<0,0>(gd,vd);
+
+ momenta_grid_domain<decltype(gd),0>(gd,mg);
+ momenta_vector<decltype(vd),0>(vd,mv);
+
+ BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
+ BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+
+ momenta_grid_domain<decltype(gd),1>(gd,mg);
+ momenta_vector<decltype(vd),1>(vd,mv);
+
+ BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
+ BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+
+ momenta_grid_domain<decltype(gd),2>(gd,mg);
+ momenta_vector<decltype(vd),2>(vd,mv);
+
+ BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
+ BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
+
+ }
+}
+
+BOOST_AUTO_TEST_CASE( int_kernel_test )
+{
+ mp4_kernel<float> mp4;
+
+ float tot = 0.0;
+
+ // Check momenta 0
+
+ tot += mp4.value(-1.3f,0);
+ tot += mp4.value(-0.3f,1);
+ tot += mp4.value(0.7f,2);
+ tot += mp4.value(1.7f,3);
+
+ BOOST_REQUIRE_CLOSE(tot,1.0f,0.001);
+
+ // Check momenta 1
+
+ tot = 0.0;
+
+ tot += -1.3f*mp4.value(-1.3f,0);
+ tot += -0.3f*mp4.value(-0.3f,1);
+ tot += 0.7f*mp4.value(0.7f,2);
+ tot += 1.7f*mp4.value(1.7f,3);
+
+ BOOST_REQUIRE_SMALL(tot,0.001f);
+
+ // Check momenta 2
+
+ tot = 0.0;
+
+ tot += (1.3f)*(1.3f)*mp4.value(-1.3f,0);
+ tot += (0.3f)*(0.3f)*mp4.value(-0.3f,1);
+ tot += (0.7f)*(0.7f)*mp4.value(0.7f,2);
+ tot += (1.7f)*(1.7f)*mp4.value(1.7f,3);
+
+ BOOST_REQUIRE_SMALL(tot,0.001f);
+
+
+ //////// Check zeta 1
+
+ tot = 0.0;
+
+ z_kernel<float,1> zk1;
+
+ tot += zk1.value(-0.3f,0);
+ tot += zk1.value(0.7f,1);
+
+ BOOST_REQUIRE_CLOSE(tot,1.0f,0.001);
+
+ //////// zeta 2 is equivalent to mp4 we do not test
+
+ //////// zeta 3
+
+ z_kernel<float,3> zk3;
+
+ tot = 0.0;
+
+ // Check momenta 0
+
+ tot += zk3.value(-2.3f,0);
+ tot += zk3.value(-1.3f,1);
+ tot += zk3.value(-0.3f,2);
+ tot += zk3.value(0.7f,3);
+ tot += zk3.value(1.7f,4);
+ tot += zk3.value(2.7f,5);
+
+ BOOST_REQUIRE_CLOSE(tot,1.0f,0.001);
+
+ // Check momenta 1
+
+ tot = 0.0;
+
+ tot += -2.3*zk3.value(-2.3f,0);
+ tot += -1.3*zk3.value(-1.3f,1);
+ tot += -0.3*zk3.value(-0.3f,2);
+ tot += 0.7*zk3.value(0.7f,3);
+ tot += 1.7*zk3.value(1.7f,4);
+ tot += 2.7*zk3.value(2.7f,5);
+
+ BOOST_REQUIRE_SMALL(tot,0.001f);
+
+ // Check momenta 2
+
+ tot = 0.0;
+
+ tot += 2.3*2.3*zk3.value(-2.3f,0);
+ tot += 1.3*1.3*zk3.value(-1.3f,1);
+ tot += 0.3*0.3*zk3.value(-0.3f,2);
+ tot += 0.7*0.7*zk3.value(0.7f,3);
+ tot += 1.7*1.7*zk3.value(1.7f,4);
+ tot += 2.7*2.7*zk3.value(2.7f,5);
+
+ BOOST_REQUIRE_SMALL(tot,0.001f);
+
+ // Check momenta 3
+
+ tot = 0.0;
+
+ tot += -2.3*-2.3*-2.3*zk3.value(-2.3f,0);
+ tot += -1.3*-1.3*-1.3*zk3.value(-1.3f,1);
+ tot += -0.3*-0.3*-0.3*zk3.value(-0.3f,2);
+ tot += 0.7*0.7*0.7*zk3.value(0.7f,3);
+ tot += 1.7*1.7*1.7*zk3.value(1.7f,4);
+ tot += 2.7*2.7*2.7*zk3.value(2.7f,5);
+
+ BOOST_REQUIRE_SMALL(tot,0.001f);
+
+
+ // z4
+
+ z_kernel<float,4> zk4;
+
+ // Check momenta 0
+
+ tot = 0.0;
+
+ tot += zk4.value(-3.3f,0);
+ tot += zk4.value(-2.3f,1);
+ tot += zk4.value(-1.3f,2);
+ tot += zk4.value(-0.3f,3);
+ tot += zk4.value(0.7f,4);
+ tot += zk4.value(1.7f,5);
+ tot += zk4.value(2.7f,6);
+ tot += zk4.value(3.7f,7);
+
+ BOOST_REQUIRE_CLOSE(tot,1.0f,0.001);
+
+ // Check momenta 1
+
+ tot = 0.0;
+
+ tot += -3.3*zk4.value(-3.3f,0);
+ tot += -2.3*zk4.value(-2.3f,1);
+ tot += -1.3*zk4.value(-1.3f,2);
+ tot += -0.3*zk4.value(-0.3f,3);
+ tot += 0.7*zk4.value(0.7f,4);
+ tot += 1.7*zk4.value(1.7f,5);
+ tot += 2.7*zk4.value(2.7f,6);
+ tot += 3.7*zk4.value(3.7f,7);
+
+ BOOST_REQUIRE_SMALL(tot,0.001f);
+
+ // Check momenta 2
+
+ tot = 0.0;
+
+ tot += 3.3*3.3*zk4.value(-3.3f,0);
+ tot += 2.3*2.3*zk4.value(-2.3f,1);
+ tot += 1.3*1.3*zk4.value(-1.3f,2);
+ tot += 0.3*0.3*zk4.value(-0.3f,3);
+ tot += 0.7*0.7*zk4.value(0.7f,4);
+ tot += 1.7*1.7*zk4.value(1.7f,5);
+ tot += 2.7*2.7*zk4.value(2.7f,6);
+ tot += 3.7*3.7*zk4.value(3.7f,7);
+
+ BOOST_REQUIRE_SMALL(tot,0.001f);
+
+ // Check momenta 3
+
+ tot = 0.0;
+
+ tot += -3.3*-3.3*-3.3*zk4.value(-3.3f,0);
+ tot += -2.3*-2.3*-2.3*zk4.value(-2.3f,1);
+ tot += -1.3*-1.3*-1.3*zk4.value(-1.3f,2);
+ tot += -0.3*-0.3*-0.3*zk4.value(-0.3f,3);
+ tot += 0.7*0.7*0.7*zk4.value(0.7f,4);
+ tot += 1.7*1.7*1.7*zk4.value(1.7f,5);
+ tot += 2.7*2.7*2.7*zk4.value(2.7f,6);
+ tot += 3.7*3.7*3.7*zk4.value(3.7f,7);
+
+ BOOST_REQUIRE_SMALL(tot,0.001f);
+
+ // Check momenta 4
+
+ tot = 0.0;
+
+ tot += -3.3*-3.3*-3.3*-3.3*zk4.value(-3.3f,0);
+ tot += -2.3*-2.3*-2.3*-2.3*zk4.value(-2.3f,1);
+ tot += -1.3*-1.3*-1.3*-1.3*zk4.value(-1.3f,2);
+ tot += -0.3*-0.3*-0.3*-0.3*zk4.value(-0.3f,3);
+ tot += 0.7*0.7*0.7*0.7*zk4.value(0.7f,4);
+ tot += 1.7*1.7*1.7*1.7*zk4.value(1.7f,5);
+ tot += 2.7*2.7*2.7*2.7*zk4.value(2.7f,6);
+ tot += 3.7*3.7*3.7*3.7*zk4.value(3.7f,7);
+
+ BOOST_REQUIRE_SMALL(tot,0.001f);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+#endif /* OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_UNIT_TESTS_HPP_ */
diff --git a/src/interpolation/z_spline.hpp b/src/interpolation/z_spline.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6988dada28f91854c9d163919f9388b12a6ffbf6
--- /dev/null
+++ b/src/interpolation/z_spline.hpp
@@ -0,0 +1,106 @@
+/*
+ * z_spline.hpp
+ *
+ * Created on: May 8, 2017
+ * Author: i-bird
+ */
+
+#ifndef OPENFPM_NUMERICS_SRC_INTERPOLATION_Z_SPLINE_HPP_
+#define OPENFPM_NUMERICS_SRC_INTERPOLATION_Z_SPLINE_HPP_
+
+#include "mp4_kernel.hpp"
+
+template<typename st, unsigned int ord>
+class z_kernel
+{
+public:
+
+ static const int np = 4;
+
+ static inline st value(st x, size_t i)
+ {
+ std::cerr << __FILE__ << ":" << __LINE__ << ": Error this order has not been implemented" << std::endl;
+ return 0.0;
+ }
+};
+
+template<typename st>
+class z_kernel<st,1>
+{
+public:
+
+ static const int np = 2;
+
+ static inline st value(st x, size_t i)
+ {
+ if (i == 0)
+ return x + 1;
+ else if (i == 1)
+ return 1-x;
+ return 0.0;
+ }
+};
+
+template<typename st>
+class z_kernel<st,2>: public mp4_kernel<st>
+{
+
+};
+
+template<typename st>
+class z_kernel<st,3>
+{
+public:
+
+ static const int np = 6;
+
+ static inline st value(st x, size_t i)
+ {
+ if (i == 0)
+ return 18.0 + ( 38.25 + (31.875 + ( 313.0/24.0 + (2.625 + 5.0/24.0*x)*x)*x)*x)*x;
+ else if (i == 1)
+ return -4.0 + (-18.75 + (-30.625 + (-545.0/24.0 + ( -7.875 - 25.0/24.0*x)*x)*x)*x)*x;
+ else if (i == 2)
+ return 1.0 + (-1.25 +(35.0/12.0 +(5.25 + 25.0/12.0*x)*x)*x)*x*x;
+ else if (i == 3)
+ return 1.0 + (-1.25 +(-35.0/12.0 +(5.25 - 25.0/12.0*x)*x)*x)*x*x;
+ else if (i == 4)
+ return -4.0 + (18.75 + (-30.625 + (545.0/24.0 + ( -7.875 + 25.0/24.0*x)*x)*x)*x)*x;
+ else if (i == 5)
+ return 18.0 + (-38.25 + (31.875 + (-313.0/24.0 + (2.625 - 5.0/24.0*x)*x)*x)*x)*x;
+
+ return 0.0;
+ }
+};
+
+template<typename st>
+class z_kernel<st,4>
+{
+public:
+
+ static const int np = 8;
+
+ static inline st value(st x, size_t i)
+ {
+ if (i == 0)
+ return 726.4 + ( 1491.2 + (58786.0/45.0 + ( 633.0 + (26383.0/144.0 + (22807.0/720.0 + (727.0/240.0 + 89.0/720.0*x)*x)*x)*x)*x)*x)*x;
+ else if (i == 1)
+ return -440 +( -1297.45 + ( -117131/72.0 + ( -1123.5 + ( -66437.0/144.0 + ( -81109.0/720.0 + ( -727.0/48.0 - 623.0/720.0*x)*x)*x)*x)*x)*x)*x;
+ else if (i == 2)
+ return 27.6 + (8617.0/60.0 +(321.825 +(395.5 +( 284.8125 + (119.7875 +(27.2625 + 623.0/240.0*x)*x)*x)*x)*x)*x)*x;
+ else if (i == 3)
+ return 1.0 + (( -49.0/36.0 + (( -959.0/144.0 + ( -2569.0/144.0 + ( -727.0/48.0 - 623.0/144.0*x)*x)*x)*x)*x)*x)*x;
+ else if (i == 4)
+ return 1.0 + (( -49.0/36.0 + (( -959.0/144.0 + ( 2569.0/144.0 + ( -727.0/48.0 + 623.0/144.0*x)*x)*x)*x)*x)*x)*x;
+ else if (i == 5)
+ return 27.6 + (-8617.0/60.0 +(321.825 +(-395.5 +( 284.8125 + (-119.7875 +(27.2625 - 623.0/240.0*x)*x)*x)*x)*x)*x)*x;
+ else if (i == 6)
+ return -440 +( 1297.45 + ( -117131/72.0 + ( 1123.5 + ( -66437.0/144.0 + ( 81109.0/720.0 + ( -727.0/48.0 + 623.0/720.0*x)*x)*x)*x)*x)*x)*x;
+ else if (i == 7)
+ return 726.4 + ( -1491.2 + (58786.0/45.0 + ( -633.0 + (26383.0/144.0 + (-22807.0/720.0 + (727.0/240.0 - 89.0/720.0*x)*x)*x)*x)*x)*x)*x;
+
+ return 0.0;
+ }
+};
+
+#endif /* OPENFPM_NUMERICS_SRC_INTERPOLATION_Z_SPLINE_HPP_ */