diff --git a/src/dinterpolate.tcc b/src/dinterpolate.tcc
index d545788..ecd04ae 100644
--- a/src/dinterpolate.tcc
+++ b/src/dinterpolate.tcc
@@ -98,7 +98,8 @@ namespace CosmoTool {
 	int signum;
 
 	gsl_linalg_LU_decomp(&M.matrix, p, &signum);
-	if (fabs(gsl_linalg_LU_det(&M.matrix, signum)) < 1e-10)
+        double a = fabs(gsl_linalg_LU_det(&M.matrix, signum));
+	if (a < 1e-10)
 	  throw InvalidArgumentException("Invalid tesselation. One simplex is coplanar.");
 	gsl_linalg_LU_invert(&M.matrix, p, &iM.matrix);
 
diff --git a/src/sphSmooth.hpp b/src/sphSmooth.hpp
index bbcdcea..b90758d 100644
--- a/src/sphSmooth.hpp
+++ b/src/sphSmooth.hpp
@@ -25,7 +25,7 @@ namespace CosmoTool
     typedef void (*runParticleValue)(ValType& t);
 
   public:
-    SPHSmooth(SPHTree *tree, uint32_t Nsph, uint32_t deltaN);
+    SPHSmooth(SPHTree *tree, uint32_t Nsph);
     virtual ~SPHSmooth();   
 
     void fetchNeighbours(const typename SPHTree::coords& c);
@@ -62,6 +62,18 @@ namespace CosmoTool
     virtual ComputePrecision getKernel(ComputePrecision d) const;    
 
     SPHTree *getTree() { return tree; }
+
+    void changeNgb(uint32_t newMax) {
+      delete[] ngb;
+      delete[] distances;
+      ngb = new SPHCell *[newMax];
+      distances = new CoordType[newMax];
+      maxNgb = newMax;
+    }
+
+    uint32_t getCurrent() const { return currentNgb; }
+
+    uint32_t getNgb() const { return maxNgb; }
  
   protected:
     SPHCell **ngb;
diff --git a/src/sphSmooth.tcc b/src/sphSmooth.tcc
index fd2f368..8d317c9 100644
--- a/src/sphSmooth.tcc
+++ b/src/sphSmooth.tcc
@@ -3,14 +3,13 @@
 namespace CosmoTool {
 
 template<typename ValType, int Ndims>
-SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph, uint32_t deltaN)
+SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph)
 {
   this->Nsph = Nsph;
-  this->deltaNsph = deltaN;
   this->tree = tree;
   this->currentNgb = 0;
 
-  maxNgb = Nsph;
+  this->maxNgb = Nsph;
   ngb = new SPHCell *[maxNgb];  
   distances = new CoordType[maxNgb];
 }
@@ -43,7 +42,6 @@ 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;
 
@@ -53,6 +51,7 @@ SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c)
   currentNgb = 0;
   for (uint32_t i = 0; i < maxNgb && ngb[i] != 0; i++,currentNgb++)
     {
+      distances[i] = sqrt(distances[i]);
       d2 = distances[i];
       if (d2 > max_dist)
 	max_dist = d2;
@@ -66,18 +65,17 @@ 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,
+  currentNgb = tree->getIntersection(c, radius, ngb, distances,
 				  maxNgb);
 
-  currentNgb = 0;
-  for (uint32_t i = 0; i < maxNgb && ngb[i] != 0; i++,currentNgb++)
+  for (uint32_t i = 0; i < currentNgb; i++)
     {
+      distances[i] = sqrt(distances[i]);
       d2 = distances[i];
       if (d2 > max_dist)
 	max_dist = d2;
@@ -119,8 +117,8 @@ ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typena
 
   for (uint32_t i = 0; i < currentNgb; i++)
     {
-      outputValue += computeWValue(c, ngb[i], distances[i], fun);
-      weight += computeWValue(c, ngb[i], distances[i], interpolateOne);
+      outputValue += computeWValue(c, *ngb[i], distances[i], fun);
+      weight += computeWValue(c, *ngb[i], distances[i], interpolateOne);
     }
 
   return (outputValue == 0) ? 0 : (outputValue / weight);