Added support for gadget2 snap

src/fortran.cpp

@@ -71,6 +71,16 @@ UnformattedRead::~UnformattedRead()
   delete f;
 }
 
+int64_t UnformattedRead::position() const
+{
+ return f->tellg();
+}
+
+void UnformattedRead::seek(int64_t pos)
+{
+  f->seekg(pos, istream::beg);
+}
+
 // Todo implement primitive description
 void UnformattedRead::setOrdering(Ordering o)
 {

src/fortran.hpp

@@ -95,6 +95,11 @@ namespace CosmoTool
     void skip(int64_t off)
       throw (InvalidUnformattedAccess);
 
+    int64_t position() const;
+    void seek(int64_t pos);
+
+    void readOrderedBuffer(void *buffer, int size)
+      throw (InvalidUnformattedAccess);
   protected:
     bool swapOrdering;
     CheckpointSize cSize;
@@ -102,8 +107,6 @@ namespace CosmoTool
     uint64_t checkPointAccum;
     std::ifstream *f;
 
-    void readOrderedBuffer(void *buffer, int size)
-      throw (InvalidUnformattedAccess);
   };
 
  class UnformattedWrite: public FortranTypes

src/loadGadget.cpp

@@ -105,6 +105,11 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
   boost::function0<float> readToSingle;
   long float_size;
 
+  if (GadgetFormat > 2) {
+     cerr << "Unknown gadget format" << endl;
+     return 0;
+  }
+
   if (id >= 0) {
     int k = snprintf(0, 0, "%s.%d", fname, id)+1;
     char *out_fname = new char[k];
@@ -124,6 +129,29 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
 
   }
 
+  typedef std::map<std::string, int64_t> BlockMap;
+  BlockMap blockTable;
+  
+  if (GadgetFormat == 2) {
+    int64_t startBlock = 0;
+    char block[5];
+    int32_t blocksize;
+
+    try {
+      while (true) {
+        f->beginCheckpoint();
+        f->readOrderedBuffer(block, 4);
+        block[4] = 0;
+        blocksize = f->readInt32();
+        f->endCheckpoint();
+        blockTable[block] = f->position();
+        f->skip(blocksize);
+      }
+    } catch (EndOfFileException&) {}
+
+    f->seek(startBlock);
+  }
+
   data = new SimuData;
   if (data == 0) {
     delete f;
@@ -131,19 +159,26 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
   }
 
   ssize_t NumPart = 0, NumPartTotal = 0;
+#define ENSURE(name) { \
+    if (GadgetFormat == 2) { \
+     BlockMap::iterator iter = blockTable.find(name); \
+     if (iter == blockTable.end()) { \
+        std::cerr << "GADGET2: Cannot find block named '" << name << "'" << endl; \
+        if (data) delete data; \
+        delete f; \
+        return 0; \
+     } \
+     f->seek(iter->second); \
+    } \
+  } 
+
   
   try
     {
+      ENSURE("HEAD");
       loadGadgetHeader(f, h, data, id);
 
-      if (GadgetFormat == 1)
       velmul = sqrt(h.time);
-      else if (GadgetFormat == 2)
-        velmul = 1/(h.time);
-      else {
-        cerr << "unknown gadget format" << endl;
-        abort();
-      }
       
       if (h.flag_doubleprecision) {
         //cout << "Gadget format with double precision" << endl;
@@ -167,6 +202,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
       return 0;
     }
 
+
   if (loadflags & NEED_TYPE)
     {
       int p = 0;
@@ -189,6 +225,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
     
     try
       {
+        ENSURE("POS ");
 	f->beginCheckpoint();
         if (f->getBlockSize() != NumPart*float_size*3) {
           // Check that single would work
@@ -239,6 +276,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
     
     try
       {
+        ENSURE("VEL ");
 	f->beginCheckpoint();
 	for(int k = 0, p = 0; k < 6; k++) {
 	  for(int n = 0; n < h.npart[k]; n++) {
@@ -271,6 +309,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
   if (loadflags & NEED_GADGET_ID) {
     try
       {
+        ENSURE("ID  ");
 	f->beginCheckpoint();
 	data->Id = new long[data->NumPart];
 	if (data->Id == 0)
@@ -314,8 +353,10 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
     try
       {
         long l = 0;
-	if (do_load)
+	if (do_load) {
+          ENSURE("MASS");
           f->beginCheckpoint();
+        }
         data->Mass = new float[NumPart];
         for (int k = 0; k < 6; k++)
 	  {
@@ -360,6 +401,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
   return data;
 }
 
+#undef ENSURE
 
 
 void CosmoTool::writeGadget(const std::string& fname, SimuData *data, int GadgetFormat)
@@ -419,7 +461,6 @@ void CosmoTool::writeGadget(const std::string& fname, SimuData *data, int Gadget
   f->endCheckpoint();
 
   float velmul = 1.0;
-  if (GadgetFormat == 1)
   velmul = sqrt(data->time);
 
   f->beginCheckpoint();