diff --git a/src/mykdtree.hpp b/src/mykdtree.hpp
index 56a9898..e87b376 100644
--- a/src/mykdtree.hpp
+++ b/src/mykdtree.hpp
@@ -8,19 +8,19 @@
 namespace CosmoTool {
 
   template<int N> 
-  struct NGBDef
+  struct KDDef
   {
     typedef float CoordType;
-    typedef float NGBCoordinates[N];
+    typedef float KDCoordinates[N];
     static const int NumCubes = 1 << (N);
   };
   
   template<int N, typename ValType>
-  struct NGBCell
+  struct KDCell
   {
     bool active;
     ValType val;
-    typename NGBDef<N>::NGBCoordinates coord;
+    typename KDDef<N>::KDCoordinates coord;
   };
 
   class NotEnoughCells: public Exception
@@ -31,11 +31,11 @@ namespace CosmoTool {
   };
 
   template<int N, typename ValType>
-  struct NGBTreeNode
+  struct KDTreeNode
   {
-    NGBCell<N,ValType> *value;
-    NGBTreeNode<N,ValType> *children[2];
-    typename NGBDef<N>::NGBCoordinates minBound, maxBound;
+    KDCell<N,ValType> *value;
+    KDTreeNode<N,ValType> *children[2];
+    typename KDDef<N>::KDCoordinates minBound, maxBound;
   };
 
   template<int N, typename ValType>
@@ -43,9 +43,10 @@ namespace CosmoTool {
   {
   public:
 
-    typename NGBDef<N>::NGBCoordinates x;
-    typename NGBDef<N>::CoordType r, r2;
-    NGBCell<N, ValType> **cells;
+    typename KDDef<N>::KDCoordinates x;
+    typename KDDef<N>::CoordType r, r2;
+    KDCell<N, ValType> **cells;
+    typename KDDef<N>::CoordType *distances;
     uint32_t currentRank;
     uint32_t numCells;
   };
@@ -54,12 +55,12 @@ namespace CosmoTool {
   class RecursionMultipleInfo
   {
   public:
-    const typename NGBDef<N>::NGBCoordinates& x;
-    BoundedQueue< NGBCell<N,ValType> *, float> queue;
+    const typename KDDef<N>::KDCoordinates& x;
+    BoundedQueue< KDCell<N,ValType> *, typename KDDef<N>::CoordType> queue;
     int traversed;
 
-    RecursionMultipleInfo(const typename NGBDef<N>::NGBCoordinates& rx,
-			  NGBCell<N,ValType> **cells,
+    RecursionMultipleInfo(const typename KDDef<N>::KDCoordinates& rx,
+			  KDCell<N,ValType> **cells,
 			  uint32_t numCells)
       : x(rx), queue(cells, numCells, INFINITY),traversed(0)
     {      
@@ -67,67 +68,76 @@ namespace CosmoTool {
   };
 
   template<int N, typename ValType>
-  class NGBTree
+  class KDTree
   {
   public:
-    static const int NumCubes = NGBDef<N>::NumCubes;
+    static const int NumCubes = KDDef<N>::NumCubes;
 
   public:
-    typedef typename NGBDef<N>::CoordType CoordType;
-    typedef typename NGBDef<N>::NGBCoordinates coords;
+    typedef typename KDDef<N>::CoordType CoordType;
+    typedef typename KDDef<N>::KDCoordinates coords;
 
-    NGBTree(NGBCell<N,ValType> *cells, uint32_t Ncells);
-    ~NGBTree();
+    KDTree(KDCell<N,ValType> *cells, uint32_t Ncells);
+    ~KDTree();
 
     uint32_t getIntersection(const coords& x, CoordType r, 
-			     NGBCell<N, ValType> **cells,
+			     KDCell<N, ValType> **cells,
+			     uint32_t numCells)
+      throw (NotEnoughCells);
+    uint32_t getIntersection(const coords& x, CoordType r, 
+			     KDCell<N, ValType> **cells,
+			     CoordType *distances,
 			     uint32_t numCells)
       throw (NotEnoughCells);
 
-    NGBCell<N, ValType> *getNearestNeighbour(const coords& x);
-    void getNearestNeighbours(const coords& x, uint32_t N,
-			      NGBCell<N, ValType> **cells);
+    KDCell<N, ValType> *getNearestNeighbour(const coords& x);
 
-    NGBTreeNode<N,ValType> *getRoot() { return root; }
+    void getNearestNeighbours(const coords& x, uint32_t N,
+			      KDCell<N, ValType> **cells);
+    void getNearestNeighbours(const coords& x, uint32_t N,
+			      KDCell<N, ValType> **cells,
+			      CoordType *distances);
+
+    KDTreeNode<N,ValType> *getRoot() { return root; }
 
     void optimize();
 
-    NGBTreeNode<N,ValType> *getAllNodes() { return nodes; }
+    KDTreeNode<N,ValType> *getAllNodes() { return nodes; }
     uint32_t getNumNodes() const { return lastNode; }
 
     uint32_t countActives() const;
 
   protected:
-    NGBTreeNode<N, ValType> *nodes;
+    KDTreeNode<N, ValType> *nodes;
     uint32_t numNodes;
     uint32_t lastNode;
 
-    NGBTreeNode<N, ValType> *root;
-    NGBCell<N, ValType> **sortingHelper;
+    KDTreeNode<N, ValType> *root;
+    KDCell<N, ValType> **sortingHelper;
 
-    NGBTreeNode<N, ValType> *buildTree(NGBCell<N,ValType> **cell0,
+    KDTreeNode<N, ValType> *buildTree(KDCell<N,ValType> **cell0,
 				       uint32_t N,
 				       uint32_t depth,
 				       coords minBound,
 				       coords maxBound);
     
     void recursiveIntersectionCells(RecursionInfoCells<N,ValType>& info,
-				    NGBTreeNode<N,ValType> *node,
+				    KDTreeNode<N,ValType> *node,
 				    int level)
       throw (NotEnoughCells);
 
-    CoordType computeDistance(NGBCell<N,ValType> *cell, const coords& x);
-    void recursiveNearest(NGBTreeNode<N, ValType> *node,
+    CoordType computeDistance(KDCell<N,ValType> *cell, const coords& x);
+    void recursiveNearest(KDTreeNode<N, ValType> *node,
 			  int level,
 			  const coords& x,
 			  CoordType& R2,
-			  NGBCell<N,ValType>*& cell);
-    void recursiveMultipleNearest(RecursionMultipleInfo<N,ValType>& info, NGBTreeNode<N,ValType> *node,
+			  KDCell<N,ValType>*& cell);
+    void recursiveMultipleNearest(RecursionMultipleInfo<N,ValType>& info, KDTreeNode<N,ValType> *node,
 				  int level);    
   };
 
   template<int N, class T>
-  uint32_t gatherActiveCells(NGBCell<N,T> **cells, uint32_t numCells);
+  uint32_t gatherActiveCells(KDCell<N,T> **cells, uint32_t numCells);
 
 };
 
diff --git a/src/mykdtree.tcc b/src/mykdtree.tcc
index 924e933..66ec29e 100644
--- a/src/mykdtree.tcc
+++ b/src/mykdtree.tcc
@@ -14,7 +14,7 @@ namespace CosmoTool {
       rank = k;
     }
 
-    bool operator()(const NGBCell<N,ValType> *a, const NGBCell<N,ValType> *b) const
+    bool operator()(const KDCell<N,ValType> *a, const KDCell<N,ValType> *b) const
     {
       return (a->coord[rank] < b->coord[rank]);
     }
@@ -23,18 +23,18 @@ namespace CosmoTool {
   };
 
   template<int N, typename ValType>
-  NGBTree<N,ValType>::~NGBTree()
+  KDTree<N,ValType>::~KDTree()
   {
   }
 
   template<int N, typename ValType>
-  NGBTree<N,ValType>::NGBTree(NGBCell<N,ValType> *cells, uint32_t Ncells)
+  KDTree<N,ValType>::KDTree(KDCell<N,ValType> *cells, uint32_t Ncells)
   {
     
     numNodes = Ncells;
-    nodes = new NGBTreeNode<N,ValType>[numNodes];
+    nodes = new KDTreeNode<N,ValType>[numNodes];
 
-    sortingHelper = new NGBCell<N,ValType> *[Ncells];
+    sortingHelper = new KDCell<N,ValType> *[Ncells];
     for (uint32_t i = 0; i < Ncells; i++)
 	sortingHelper[i] = &cells[i];    
 
@@ -42,7 +42,7 @@ namespace CosmoTool {
   }
 
   template<int N, typename ValType>
-  void NGBTree<N,ValType>::optimize()
+  void KDTree<N,ValType>::optimize()
   {
     coords absoluteMin, absoluteMax;
 
@@ -63,8 +63,8 @@ namespace CosmoTool {
   }
 
   template<int N, typename ValType>
-  uint32_t NGBTree<N,ValType>::getIntersection(const coords& x, CoordType r, 
-					       NGBCell<N, ValType> **cells,
+  uint32_t KDTree<N,ValType>::getIntersection(const coords& x, CoordType r, 
+					       KDCell<N, ValType> **cells,
 					       uint32_t numCells)
     throw (NotEnoughCells)
   {
@@ -76,14 +76,36 @@ namespace CosmoTool {
     info.cells = cells;
     info.currentRank = 0;
     info.numCells = numCells;
+    info.distances = 0;
 
     recursiveIntersectionCells(info, root, 0);
     return info.currentRank;
   }
 
   template<int N, typename ValType>
-  void NGBTree<N,ValType>::recursiveIntersectionCells(RecursionInfoCells<N,ValType>& info, 
-						      NGBTreeNode<N,ValType> *node,
+  uint32_t KDTree<N,ValType>::getIntersection(const coords& x, CoordType r, 
+					      KDCell<N, ValType> **cells,
+					      CoordType *distances,
+					      uint32_t numCells)
+    throw (NotEnoughCells)
+  {
+    RecursionInfoCells<N,ValType> info;
+
+    memcpy(info.x, x, sizeof(x));
+    info.r = r;
+    info.r2 = r*r;
+    info.cells = cells;
+    info.currentRank = 0;
+    info.numCells = numCells;
+    info.distances = distances;
+
+    recursiveIntersectionCells(info, root, 0);
+    return info.currentRank;
+  }
+
+  template<int N, typename ValType>
+  void KDTree<N,ValType>::recursiveIntersectionCells(RecursionInfoCells<N,ValType>& info, 
+						      KDTreeNode<N,ValType> *node,
 						      int level)
     throw (NotEnoughCells)
   {
@@ -101,7 +123,10 @@ namespace CosmoTool {
 	  {
 	    if (info.currentRank == info.numCells)
 	      throw NotEnoughCells();
-	    info.cells[info.currentRank++] = node->value;
+	    info.cells[info.currentRank] = node->value;
+	    if (info.distances)
+	      info.distances[info.currentRank] = d2;
+	    info.currentRank++;
 	  }
       }
 
@@ -123,7 +148,7 @@ namespace CosmoTool {
   }
 
   template<int N, typename ValType>
-  uint32_t gatherActiveCells(NGBCell<N,ValType> **cells,
+  uint32_t gatherActiveCells(KDCell<N,ValType> **cells,
 			 uint32_t Ncells)
   {
     uint32_t swapId = Ncells-1;
@@ -148,7 +173,7 @@ namespace CosmoTool {
   }
 
   template<int N, typename ValType>
-  NGBTreeNode<N, ValType> *NGBTree<N,ValType>::buildTree(NGBCell<N,ValType> **cell0,
+  KDTreeNode<N, ValType> *KDTree<N,ValType>::buildTree(KDCell<N,ValType> **cell0,
 							 uint32_t Ncells,
 							 uint32_t depth,
 							 coords minBound,
@@ -158,7 +183,7 @@ namespace CosmoTool {
       return 0;
 
     int axis = depth % N;
-    NGBTreeNode<N,ValType> *node = &nodes[lastNode++];
+    KDTreeNode<N,ValType> *node = &nodes[lastNode++];
     uint32_t mid = Ncells/2;
     coords tmpBound;
     
@@ -187,7 +212,7 @@ namespace CosmoTool {
   }
 
   template<int N, typename ValType>
-  uint32_t NGBTree<N,ValType>::countActives() const
+  uint32_t KDTree<N,ValType>::countActives() const
   {
     uint32_t numActive = 0;
     for (uint32_t i = 0; i < lastNode; i++)
@@ -199,8 +224,8 @@ namespace CosmoTool {
   }
 
   template<int N, typename ValType>
-  typename NGBDef<N>::CoordType
-  NGBTree<N,ValType>::computeDistance(NGBCell<N, ValType> *cell, const coords& x)
+  typename KDDef<N>::CoordType
+  KDTree<N,ValType>::computeDistance(KDCell<N, ValType> *cell, const coords& x)
   {
     CoordType d2 = 0;
 
@@ -214,16 +239,16 @@ namespace CosmoTool {
 
   template<int N, typename ValType>
   void
-  NGBTree<N,ValType>::recursiveNearest(
-				       NGBTreeNode<N,ValType> *node,
+  KDTree<N,ValType>::recursiveNearest(
+				       KDTreeNode<N,ValType> *node,
 				       int level,
 				       const coords& x,
 				       CoordType& R2,
-				       NGBCell<N,ValType> *& best)
+				       KDCell<N,ValType> *& best)
   {
     CoordType d2 = 0;
     int axis = level % N;
-    NGBTreeNode<N,ValType> *other, *go;
+    KDTreeNode<N,ValType> *other, *go;
 
     if (x[axis] < node->value->coord[axis])
       {
@@ -287,11 +312,11 @@ namespace CosmoTool {
   }
 
  template<int N, typename ValType>
- NGBCell<N, ValType> *
- NGBTree<N,ValType>::getNearestNeighbour(const coords& x)
+ KDCell<N, ValType> *
+ KDTree<N,ValType>::getNearestNeighbour(const coords& x)
  {
    CoordType R2 = INFINITY;
-   NGBCell<N,ValType> *best = 0;
+   KDCell<N,ValType> *best = 0;
 
    recursiveNearest(root, 0, x, R2, best);   
    
@@ -300,12 +325,12 @@ namespace CosmoTool {
 
   template<int N, typename ValType>
   void
-  NGBTree<N,ValType>::recursiveMultipleNearest(RecursionMultipleInfo<N,ValType>& info, NGBTreeNode<N,ValType> *node,
+  KDTree<N,ValType>::recursiveMultipleNearest(RecursionMultipleInfo<N,ValType>& info, KDTreeNode<N,ValType> *node,
 					       int level)
   {
     CoordType d2 = 0;
     int axis = level % N;
-    NGBTreeNode<N,ValType> *other, *go;
+    KDTreeNode<N,ValType> *other, *go;
 
     if (info.x[axis] < node->value->coord[axis])
       {
@@ -356,8 +381,8 @@ namespace CosmoTool {
   }
 
  template<int N, typename ValType>
- void NGBTree<N,ValType>::getNearestNeighbours(const coords& x, uint32_t N2,
-					       NGBCell<N, ValType> **cells)
+ void KDTree<N,ValType>::getNearestNeighbours(const coords& x, uint32_t N2,
+					       KDCell<N, ValType> **cells)
  {
    RecursionMultipleInfo<N,ValType> info(x, cells, N2);
    
@@ -369,4 +394,20 @@ namespace CosmoTool {
    std::cout << "Traversed = " << info.traversed << std::endl;
  }
 
+ template<int N, typename ValType>
+ void KDTree<N,ValType>::getNearestNeighbours(const coords& x, uint32_t N2,
+					      KDCell<N, ValType> **cells,
+					      CoordType *distances)
+ {
+   RecursionMultipleInfo<N,ValType> info(x, cells, N2);
+   
+   for (int i = 0; i < N2; i++)
+     cells[i] = 0;
+
+   recursiveMultipleNearest(info, root, 0);
+   memcpy(distances, info.getPriorities(), sizeof(CoordType)*N2);
+
+   std::cout << "Traversed = " << info.traversed << std::endl;
+ }
+
 };