diff --git a/src/sphSmooth.hpp b/src/sphSmooth.hpp
new file mode 100644
index 0000000..bbcdcea
--- /dev/null
+++ b/src/sphSmooth.hpp
@@ -0,0 +1,92 @@
+#ifndef __COSMOTOOL_SPH_SMOOTH_HPP
+#define __COSMOTOOL_SPH_SMOOTH_HPP
+
+#include "config.hpp"
+#include "mykdtree.hpp"
+
+namespace CosmoTool
+{
+  template <typename ValType, int Ndims = NUMDIMS>
+  class SPHSmooth
+  {
+  public:
+    typedef struct
+    {
+      ComputePrecision weight;
+      ValType pValue;
+    } FullType;
+
+    typedef KDTree<Ndims,FullType> SPHTree;
+    typedef KDTreeNode<Ndims,FullType> SPHNode;
+    typedef KDCell<Ndims,FullType>  SPHCell;
+    typedef typename KDTree<Ndims,FullType>::CoordType CoordType;
+
+    typedef ComputePrecision (*computeParticleValue)(const ValType& t);
+    typedef void (*runParticleValue)(ValType& t);
+
+  public:
+    SPHSmooth(SPHTree *tree, uint32_t Nsph, uint32_t deltaN);
+    virtual ~SPHSmooth();   
+
+    void fetchNeighbours(const typename SPHTree::coords& c);
+    void fetchNeighboursOnVolume(const typename SPHTree::coords& c, ComputePrecision radius);
+    const typename SPHTree::coords& getCurrentCenter() const
+    {
+      return currentCenter;
+    }
+
+    ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c,
+					  computeParticleValue fun);    
+    ComputePrecision computeInterpolatedValue(const typename SPHTree::coords& c,
+					      computeParticleValue fun);    
+    ComputePrecision getMaxDistance(const typename SPHTree::coords& c, 
+				    SPHNode *node) const;
+
+    ComputePrecision getSmoothingLen() const
+    {
+      return smoothRadius;
+    }
+
+    // TO USE WITH EXTREME CARE !
+    void setSmoothingLen(ComputePrecision len)
+    {
+      smoothRadius = len;
+    }
+    // END
+
+    void runForEachNeighbour(runParticleValue fun);
+    void addGridSite(const typename SPHTree::coords& c);
+
+    bool hasNeighbours() const;
+
+    virtual ComputePrecision getKernel(ComputePrecision d) const;    
+
+    SPHTree *getTree() { return tree; }
+ 
+  protected:
+    SPHCell **ngb;
+    CoordType *distances;
+    uint32_t Nsph;
+    uint32_t deltaNsph;
+    uint32_t maxNgb;
+    uint32_t currentNgb;
+    SPHTree *tree;
+    ComputePrecision smoothRadius;
+    typename SPHTree::coords currentCenter;
+    
+    ComputePrecision computeWValue(const typename SPHTree::coords & c,
+				   SPHCell& cell,
+				   CoordType d,
+				   computeParticleValue fun);
+    void runUnrollNode(SPHNode *node,
+		       runParticleValue fun);
+  };
+
+  template<class ValType1, class ValType2, int Ndims>
+  bool operator<(const SPHSmooth<ValType1, Ndims>& s1, const SPHSmooth<ValType2, Ndims>& s2);
+
+};
+
+#include "sphSmooth.tcc"
+
+#endif
diff --git a/src/sphSmooth.tcc b/src/sphSmooth.tcc
new file mode 100644
index 0000000..fd2f368
--- /dev/null
+++ b/src/sphSmooth.tcc
@@ -0,0 +1,184 @@
+#include <cmath>
+
+namespace CosmoTool {
+
+template<typename ValType, int Ndims>
+SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph, uint32_t deltaN)
+{
+  this->Nsph = Nsph;
+  this->deltaNsph = deltaN;
+  this->tree = tree;
+  this->currentNgb = 0;
+
+  maxNgb = Nsph;
+  ngb = new SPHCell *[maxNgb];  
+  distances = new CoordType[maxNgb];
+}
+
+template<typename ValType, int Ndims>
+SPHSmooth<ValType,Ndims>::~SPHSmooth()
+{
+  delete[] ngb;
+  delete[] distances;
+}
+
+template<typename ValType, int Ndims>
+ComputePrecision SPHSmooth<ValType,Ndims>::computeWValue(const typename SPHTree::coords& c,
+							 SPHCell& cell,
+							 CoordType d,
+							 computeParticleValue fun)
+{
+  CoordType weight;
+
+  d /= smoothRadius;
+  weight = getKernel(d);
+
+  if (cell.val.weight != 0)
+    return weight * fun(cell.val.pValue) / cell.val.weight;
+  else
+    return 0;
+}
+
+template<typename ValType, int Ndims>
+void
+SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c)
+{
+  typename SPHTree::coords radius;
+  ComputePrecision d2, max_dist = 0;
+  uint32_t requested = Nsph;
+
+  memcpy(currentCenter, c, sizeof(c));
+  tree->getNearestNeighbours(c, maxNgb, ngb, distances);
+
+  currentNgb = 0;
+  for (uint32_t i = 0; i < maxNgb && ngb[i] != 0; i++,currentNgb++)
+    {
+      d2 = distances[i];
+      if (d2 > max_dist)
+	max_dist = d2;
+    }
+
+  smoothRadius = max_dist / 2;  
+}					
+
+template<typename ValType, int Ndims>
+void
+SPHSmooth<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords& c,
+						  ComputePrecision radius)
+{
+  typename SPHTree::coords allRadii = {radius,radius,radius};
+  uint32_t numPart;
+  ComputePrecision d2, max_dist = 0;
+
+  memcpy(currentCenter, c, sizeof(c));
+
+  numPart = tree->getIntersection(c, ngb, distances,
+				  maxNgb);
+
+  currentNgb = 0;
+  for (uint32_t i = 0; i < maxNgb && ngb[i] != 0; i++,currentNgb++)
+    {
+      d2 = distances[i];
+      if (d2 > max_dist)
+	max_dist = d2;
+    }
+  smoothRadius = max_dist / 2;  
+}				
+
+template<typename ValType, int Ndims>
+ComputePrecision 
+SPHSmooth<ValType,Ndims>::computeSmoothedValue(const typename SPHTree::coords& c,
+					       computeParticleValue fun)
+{
+  ComputePrecision outputValue = 0;
+  ComputePrecision max_dist = 0;
+  ComputePrecision r3 = smoothRadius * smoothRadius * smoothRadius;
+
+  for (uint32_t i = 0; i < currentNgb; i++)
+    {
+      outputValue += computeWValue(c, *ngb[i], distances[i], fun);
+    }
+
+  return outputValue / r3;
+}
+
+template<typename ValType>
+ComputePrecision interpolateOne(const ValType& t)
+{
+  return 1.0;
+}
+
+// WARNING ! Cell's weight must be 1 !!!
+template<typename ValType, int Ndims>
+ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typename SPHTree::coords& c,
+								    computeParticleValue fun)
+{
+  ComputePrecision outputValue = 0;
+  ComputePrecision max_dist = 0;
+  ComputePrecision weight = 0;
+
+  for (uint32_t i = 0; i < currentNgb; i++)
+    {
+      outputValue += computeWValue(c, ngb[i], distances[i], fun);
+      weight += computeWValue(c, ngb[i], distances[i], interpolateOne);
+    }
+
+  return (outputValue == 0) ? 0 : (outputValue / weight);  
+}
+
+template<typename ValType, int Ndims>
+void SPHSmooth<ValType,Ndims>::runForEachNeighbour(runParticleValue fun)
+{
+  for (uint32_t i = 0; i < currentNgb; i++)
+    {
+      fun(ngb[i]);
+    }
+}
+
+
+template<typename ValType, int Ndims>
+void SPHSmooth<ValType,Ndims>::addGridSite(const typename SPHTree::coords& c)
+{
+  ComputePrecision outputValue = 0;
+  ComputePrecision max_dist = 0;
+  
+  ComputePrecision r3 = smoothRadius * smoothRadius * smoothRadius;
+
+  for (uint32_t i = 0; i < currentNgb; i++)
+    {
+      ComputePrecision d = distances[i];
+      SPHCell& cell = *ngb[i];
+      cell.val.weight += getKernel(d/smoothRadius) / r3;
+
+    }
+}
+
+template<typename ValType, int Ndims>
+ComputePrecision 
+SPHSmooth<ValType,Ndims>::getKernel(ComputePrecision x) const
+{
+  // WARNING !!! This is an unnormalized version of the kernel.
+  if (x < 1)
+    return 1 - 1.5 * x * x + 0.75 * x * x * x;
+  else if (x < 2)
+    {
+      ComputePrecision d = 2 - x;
+      return 0.25 * d * d * d;
+    }
+  else
+    return 0;
+}
+
+template<typename ValType, int Ndims>
+bool SPHSmooth<ValType,Ndims>::hasNeighbours() const
+{
+  return (currentNgb != 0);
+}
+
+template<class ValType1, class ValType2, int Ndims>
+bool operator<(const SPHSmooth<ValType1, Ndims>& s1, const SPHSmooth<ValType2, Ndims>& s2)
+{
+  return (s1.getSmoothingLen() < s2.getSmoothingLen());
+}
+
+};