Fixed to correctly support velocities in gadget. Fixed on-disk format on kdtree

src/loadGadget.cpp

@@ -1,3 +1,4 @@
+#include <cmath>
 #include <iostream>
 #include <assert.h>
 #include <stdio.h>
@@ -63,12 +64,13 @@ PurePositionData *CosmoTool::loadGadgetPosition(const char *fname)
 
 
 
-SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags)
+SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, int GadgetFormat)
 {
   SimuData *data;
   int p, n;
   UnformattedRead *f;
   GadgetHeader h;
+  float velmul;
 
   if (id >= 0) {
     int k = snprintf(0, 0, "%s.%d", fname, id)+1;
@@ -113,9 +115,9 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags)
       h.flag_cooling = f->readInt32();
       h.num_files = f->readInt32();
       data->BoxSize = h.BoxSize = f->readReal64();
-      h.Omega0 = f->readReal64();
+      data->Omega_M = h.Omega0 = f->readReal64();
-      h.OmegaLambda = f->readReal64();
+      data->Omega_Lambda = h.OmegaLambda = f->readReal64();
-      h.HubbleParam = f->readReal64();
+      data->Hubble = h.HubbleParam = f->readReal64();
       f->endCheckpoint(true);
       
       for(int k=0; k<6; k++)
@@ -125,6 +127,14 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags)
 	}
       data->NumPart = NumPart;
       data->TotalNumPart = NumPartTotal;
+      if (GadgetFormat == 1)
+        velmul = sqrt(h.time);
+      else if (GadgetFormat == 2)
+        velmul = 1/(h.time);
+      else {
+        cerr << "unknown gadget format" << endl;
+	abort();
+      }
     }
   catch (const InvalidUnformattedAccess& e)
     {
@@ -197,9 +207,10 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags)
 	f->beginCheckpoint();
 	for(int k = 0, p = 0; k < 6; k++) {
 	  for(int n = 0; n < h.npart[k]; n++) {
-	    data->Vel[0][p] = f->readReal32();
+	    // THIS IS GADGET 1
-	    data->Vel[1][p] = f->readReal32();
+	    data->Vel[0][p] = f->readReal32()*velmul;
-	    data->Vel[2][p] = f->readReal32();
+	    data->Vel[1][p] = f->readReal32()*velmul;
+	    data->Vel[2][p] = f->readReal32()*velmul;
 	    p++;
 	  }
 	}
@@ -213,7 +224,8 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags)
 	return 0;
       }
     
-    // TODO: FIX THE UNITS OF THESE FUNKY VELOCITIES !!!
+    // THE VELOCITIES ARE IN PHYSICAL COORDINATES
+///    // TODO: FIX THE UNITS OF THESE FUNKY VELOCITIES !!!
   } else {
     // Skip velocities
     f->skip(NumPart*3*sizeof(float)+2*4);

src/loadGadget.hpp

@@ -8,7 +8,7 @@ namespace CosmoTool {
   
   PurePositionData *loadGadgetPosition(const char *fname);
 
-  SimuData *loadGadgetMulti(const char *fname, int id, int flags);
+  SimuData *loadGadgetMulti(const char *fname, int id, int flags, int GadgetFormat = 1);
   
 };
 

src/loadRamses.cpp

@@ -243,7 +243,6 @@ int readInfoFile(const char *basename, int outputId, InfoData& info)
 	  iss >> unitLength;
 
 	  infoMap[keyword] = unitLength;
-	  
 	}
     }
   

src/mykdtree.tcc

@@ -446,6 +446,7 @@ namespace CosmoTool {
   {
     char id[KDTREE_DISK_SIGNATURE_LEN];
     long nodesUsed, numCells;
+    long rootId;
   };
 
   template<int N, typename ValType, typename CType>
@@ -456,6 +457,7 @@ namespace CosmoTool {
     strncpy(h.id, KDTREE_DISK_SIGNATURE, KDTREE_DISK_SIGNATURE_LEN);
     h.nodesUsed = lastNode;
     h.numCells = numNodes;
+    h.rootId = root - nodes;
     o.write((char*)&h, sizeof(h));
 
     for (long i = 0; i < lastNode; i++)
@@ -464,13 +466,17 @@ namespace CosmoTool {
 
 	node_on_disk.cell_id = nodes[i].value - base_cell;
 	if (nodes[i].children[0] == 0)
-	  node_on_disk.children_node[0] = -1;
+	  node_on_disk.children_node[0] = -1L;
 	else
 	  node_on_disk.children_node[0] = nodes[i].children[0] - nodes;
+       assert((node_on_disk.children_node[0] == -1) || ((node_on_disk.children_node[0] >= 0) &&  (node_on_disk.children_node[0] < lastNode)));
+
 	if (nodes[i].children[1] == 0)
-	  node_on_disk.children_node[1] = -1;
+	  node_on_disk.children_node[1] = -1L;
 	else
 	  node_on_disk.children_node[1] = nodes[i].children[1] - nodes;
+       assert((node_on_disk.children_node[1] == -1) || ((node_on_disk.children_node[1] >= 0) &&  (node_on_disk.children_node[1] < lastNode)));
+
 	memcpy(node_on_disk.minBound, nodes[i].minBound, sizeof(coords));
 	memcpy(node_on_disk.maxBound, nodes[i].maxBound, sizeof(coords));
 
@@ -489,10 +495,15 @@ namespace CosmoTool {
 
     in.read((char *)&h, sizeof(h));
     if (!in || strncmp(h.id, KDTREE_DISK_SIGNATURE, KDTREE_DISK_SIGNATURE_LEN) != 0)
-      throw InvalidOnDiskKDTree();
+      {
+         std::cerr << "KDTree Signature invalid" << std::endl;
+         throw InvalidOnDiskKDTree();
+      } 
 
-    if (h.numCells != Ncells || h.nodesUsed < 0)
+    if (h.numCells != Ncells || h.nodesUsed < 0) {
+      std::cerr << "The number of cells has changed (" << h.numCells << " != " << Ncells << ") or nodesUsed=" << h.nodesUsed << std::endl;
       throw InvalidOnDiskKDTree();
+    }
 
     base_cell = cells;
     nodes = new Node[h.nodesUsed];
@@ -505,11 +516,18 @@ namespace CosmoTool {
 
 	in.read((char *)&node_on_disk, sizeof(node_on_disk));
 
-	if (!in || node_on_disk.cell_id > numNodes || node_on_disk.cell_id < 0 ||
+	if (!in) {
+           std::cerr << "End-of-file reached" << std::endl;
+           delete[] nodes;
+           throw InvalidOnDiskKDTree();
+        }
+        if (node_on_disk.cell_id > numNodes || node_on_disk.cell_id < 0 ||
 	    node_on_disk.children_node[0] > lastNode || node_on_disk.children_node[0] < -1 ||
 	    node_on_disk.children_node[1] > lastNode || node_on_disk.children_node[1] < -1)
 	  {
 	    delete[] nodes;
+            std::cerr << "Invalid cell id or children node id invalid" << std::endl;
+            std::cerr << node_on_disk.cell_id << std::endl << node_on_disk.children_node[0] << std::endl << node_on_disk.children_node[1] << std::endl;
 	    throw InvalidOnDiskKDTree();
 	  }
 
@@ -535,10 +553,16 @@ namespace CosmoTool {
 	if (c != N)
 	  {
 	    delete[] nodes;
+            std::cerr << "Coordinates of the cell inconsistent with the boundaries" << std::endl
+                      << "   X=" << nodes[i].value->coord[0] << " B=[" << nodes[i].minBound[0] << "," << nodes[i].maxBound[0] << "]" << std::endl
+                      << "   Y=" << nodes[i].value->coord[1] << " B=[" << nodes[i].minBound[1] << "," << nodes[i].maxBound[1] << "]" << std::endl
+                      << "   Z=" << nodes[i].value->coord[2] << " B=[" << nodes[i].minBound[2] << "," << nodes[i].maxBound[2] << "]" << std::endl;
 	    throw InvalidOnDiskKDTree();
 	  }
       }
 
+    root = &nodes[h.rootId];
+    
     sortingHelper = new Cell *[Ncells];
     for (uint32_t i = 0; i < Ncells; i++)
       sortingHelper[i] = &cells[i];