From 51c5e5dcb23c1b3d97b186372ae7588018db90a0 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@curvaton.astro.uiuc.edu>
Date: Wed, 22 Sep 2010 15:33:06 -0500
Subject: [PATCH] Better ramses loading

---
 src/loadRamses.cpp | 174 ++++++++++++++++++++++++++++++++++++++++++++-
 src/loadRamses.hpp |   3 +
 2 files changed, 176 insertions(+), 1 deletion(-)

diff --git a/src/loadRamses.cpp b/src/loadRamses.cpp
index 01cd4c3..5306536 100644
--- a/src/loadRamses.cpp
+++ b/src/loadRamses.cpp
@@ -203,7 +203,7 @@ int readInfoFile(const char *basename, int outputId, InfoData& info)
   ostringstream ss_fname;  
   ss_fname << basename << "/info_" << setfill('0') << setw(5) << outputId << ".txt";
 
-  cout << "Opening info file " << ss_fname.str() << endl;
+//  cout << "Opening info file " << ss_fname.str() << endl;
   ifstream infile(ss_fname.str().c_str());
   if (!infile)
     return 0;
@@ -255,6 +255,178 @@ int readInfoFile(const char *basename, int outputId, InfoData& info)
   return 1;
 }
 
+CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputId, int cpuid, int flags)
+{
+  CosmoTool::SimuData *data = new CosmoTool::SimuData();
+
+  int id = 1;
+  uint32_t totPart = 0;
+  int nCpu = 0;
+  InfoData info;
+
+  static const double CM_IN_MPC = 3.08e24;
+
+  if (!readInfoFile(basename, outputId, info))
+    return 0;
+
+  double hubble = info.aexp*info.aexp/info.unit_t / (1e5/CM_IN_MPC);
+  double L0 = info.boxSize*info.unitLength*hubble/100/CM_IN_MPC/info.aexp;
+  double unit_vel = L0*hubble/info.aexp;
+  
+  while (1)
+    {
+      ostringstream ss_fname;
+      ss_fname << basename << "/part_" << setfill('0') << setw(5) << outputId << ".out" << setfill('0') << setw(5) << id;
+      string fname = ss_fname.str();
+
+      try
+        {
+          UnformattedRead infile(fname);
+
+          int ndim, nPar;
+
+          infile.beginCheckpoint();
+          nCpu = max(1,infile.readInt32());
+          infile.endCheckpoint();
+
+          infile.beginCheckpoint();
+          ndim = infile.readInt32();
+          infile.endCheckpoint();
+
+          infile.beginCheckpoint();
+          nPar = infile.readInt32();
+          infile.endCheckpoint();
+
+          totPart += nPar;
+          id++;
+        }
+      catch (const NoSuchFileException& e)
+        {
+          break;
+        }
+    }
+
+
+  data->BoxSize = L0;
+  data->time = info.aexp;
+  data->NumPart = 0;
+  data->TotalNumPart = totPart;
+
+  if (cpuid < 0)
+    cpuid = 1;
+  else cpuid++;
+
+  uint32_t curPos = 0;
+  
+  ostringstream ss_fname;
+  ss_fname << basename << "/part_" << setfill('0') << setw(5) << outputId << ".out" << setfill('0') << setw(5) << cpuid;
+  
+  string fname = ss_fname.str();
+    
+  try
+    {
+      UnformattedRead infile(fname);
+      
+      int nCpu, ndim, nPar;
+      
+      infile.beginCheckpoint();
+      nCpu = infile.readInt32();
+      infile.endCheckpoint();
+      
+      infile.beginCheckpoint();
+      ndim = infile.readInt32();
+      infile.endCheckpoint();
+      
+      infile.beginCheckpoint();
+      data->NumPart = nPar = infile.readInt32();
+      infile.endCheckpoint();
+      
+      if (flags & NEED_POSITION)
+	{
+	  data->Pos[0] = new float[nPar];
+	  data->Pos[1] = new float[nPar];
+	  data->Pos[2] = new float[nPar];
+	}
+      if (flags & NEED_VELOCITY)
+	{
+	  data->Vel[0] = new float[nPar];
+	  data->Vel[1] = new float[nPar];
+	  data->Vel[2] = new float[nPar];
+	}
+      if (flags & NEED_GADGET_ID)
+	{
+	  data->Id = new int[nPar];
+	}
+
+      for (int k = 0; k < 3; k++)
+	{
+	  infile.beginCheckpoint();
+	  if (flags & NEED_POSITION)
+	    {
+	      for (uint32_t i = 0; i < nPar; i++)
+		{
+		  data->Pos[k][i] = infile.readReal32();
+		  data->Pos[k][i] *= data->BoxSize;
+		}
+	      infile.endCheckpoint();
+	    }
+	  else
+	    infile.endCheckpoint(true);
+	}
+      
+      for (int k = 0; k < 3; k++) {
+	infile.beginCheckpoint();
+	if (flags & NEED_VELOCITY)
+	  {
+	    for (uint32_t i = 0; i < nPar; i++)
+	      {
+		data->Vel[k][i] = infile.readReal32();
+		data->Vel[k][i] *= unit_vel;
+	      }
+	    infile.endCheckpoint();
+	  }
+	else
+	  infile.endCheckpoint(true);
+      }
+      
+      float minMass = INFINITY;
+      infile.beginCheckpoint();
+      for (uint32_t i = nPar; i > 0; i--)
+	{
+	  float dummyF = infile.readReal32();
+	  if (dummyF < minMass) minMass = dummyF;
+	}
+      infile.endCheckpoint();
+      
+      infile.beginCheckpoint();
+      if (flags & NEED_GADGET_ID)
+	{
+	  for (uint32_t i = 0; i < nPar; i++)
+	    data->Id[i] = infile.readInt32();
+	  
+	  infile.endCheckpoint();
+	}
+      else
+	infile.endCheckpoint(true);
+      
+      curPos += nPar;
+    }
+  catch (const NoSuchFileException& e)  
+    {
+	cerr << "No such file " << fname << endl;
+      delete data;
+      return 0;
+    }
+  
+  return data;
+}
+
+
+
+
+
+
+
 CosmoTool::PurePositionData *CosmoTool::loadRamsesPosition(const char *basename, int outputId, bool quiet)
 {
   PurePositionData *gd = (PurePositionData *)malloc(sizeof(PurePositionData));
diff --git a/src/loadRamses.hpp b/src/loadRamses.hpp
index d725a70..c3efc84 100644
--- a/src/loadRamses.hpp
+++ b/src/loadRamses.hpp
@@ -2,6 +2,7 @@
 #define _LOAD_RAMSES_HPP
 
 #include "load_data.hpp"
+#include "loadSimu.hpp"
 
 namespace CosmoTool {
 
@@ -10,6 +11,8 @@ namespace CosmoTool {
   PhaseSpaceData *loadRamsesPhase(const char *fname, int id, bool quiet = false);
 
   PhaseSpaceDataID *loadRamsesPhase1(const char *fname, int id, int cpuid, bool quiet = false);
+
+  SimuData *loadRamsesSimu(const char *basename, int id, int cpuid, int flags);
 };
 
 #endif