diff --git a/src/sphSmooth.hpp b/src/sphSmooth.hpp
index b90758d..7d14b8a 100644
--- a/src/sphSmooth.hpp
+++ b/src/sphSmooth.hpp
@@ -29,6 +29,7 @@ namespace CosmoTool
     virtual ~SPHSmooth();   
 
     void fetchNeighbours(const typename SPHTree::coords& c);
+    void fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNsph);
     void fetchNeighboursOnVolume(const typename SPHTree::coords& c, ComputePrecision radius);
     const typename SPHTree::coords& getCurrentCenter() const
     {
diff --git a/src/sphSmooth.tcc b/src/sphSmooth.tcc
index 8d317c9..0cf7de8 100644
--- a/src/sphSmooth.tcc
+++ b/src/sphSmooth.tcc
@@ -38,6 +38,38 @@ ComputePrecision SPHSmooth<ValType,Ndims>::computeWValue(const typename SPHTree:
     return 0;
 }
 
+template<typename ValType, int Ndims>
+void
+SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNngb)
+{
+  ComputePrecision d2, max_dist = 0;
+  uint32_t requested = newNngb;
+
+  if (requested > maxNgb)
+    {
+      delete[] ngb;
+      delete[] distances;
+
+      maxNgb = requested;
+      ngb = new SPHCell *[maxNgb];
+      distances = new CoordType[maxNgb];
+    }
+
+  memcpy(currentCenter, c, sizeof(c));
+  tree->getNearestNeighbours(c, requested, ngb, distances);
+
+  currentNgb = 0;
+  for (uint32_t i = 0; i < requested && ngb[i] != 0; i++,currentNgb++)
+    {
+      distances[i] = sqrt(distances[i]);
+      d2 = distances[i];
+      if (d2 > max_dist)
+        max_dist = d2;
+    }
+
+  smoothRadius = max_dist / 2;  
+}
+
 template<typename ValType, int Ndims>
 void
 SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c)
@@ -46,10 +78,10 @@ SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c)
   uint32_t requested = Nsph;
 
   memcpy(currentCenter, c, sizeof(c));
-  tree->getNearestNeighbours(c, maxNgb, ngb, distances);
+  tree->getNearestNeighbours(c, requested, ngb, distances);
 
   currentNgb = 0;
-  for (uint32_t i = 0; i < maxNgb && ngb[i] != 0; i++,currentNgb++)
+  for (uint32_t i = 0; i < requested && ngb[i] != 0; i++,currentNgb++)
     {
       distances[i] = sqrt(distances[i]);
       d2 = distances[i];