Load gadget data

lib/Makefile

@@ -1,5 +1,5 @@
 SHLIBS= libCosmoTool.so
-SOURCES= loadRamses.cpp yorick.cpp miniargs.cpp  fortran.cpp interpolate.cpp load_data.cpp powerSpectrum.cpp octTree.cpp 
+SOURCES= loadRamses.cpp yorick.cpp miniargs.cpp  fortran.cpp interpolate.cpp load_data.cpp powerSpectrum.cpp octTree.cpp loadGadget.cpp
 LIBS= -lnetcdf_c++ -lnetcdf -lgsl -lgslcblas -lm
 
 include config.mk
@@ -8,7 +8,7 @@ VPATH=../src
 
 all: $(SHLIBS)
 
-libCosmoTool.so: loadRamses.o yorick.o miniargs.o fortran.o interpolate.o load_data.o powerSpectrum.o octTree.o
+libCosmoTool.so: loadRamses.o yorick.o miniargs.o fortran.o interpolate.o load_data.o powerSpectrum.o octTree.o loadGadget.o
 
 depend: $(SOURCES)
 	@echo "[DEPENDS] $^"

src/loadGadget.cpp

@@ -0,0 +1,59 @@
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "load_data.hpp"
+#include "loadGadget.hpp"
+#include "fortran.hpp"
+
+using namespace CosmoTool;
+
+PurePositionData *CosmoTool::loadGadgetPosition(const char *fname)
+{
+  PurePositionData *data;
+  int p, n;
+  UnformattedRead f(fname);
+  GadgetHeader h;
+  
+  data = new PurePositionData;
+  f.beginCheckpoint();
+  for (int i = 0; i < 6; i++)
+    h.npart[i] = f.readInt32();
+  for (int i = 0; i < 6; i++)    
+    h.mass[i] = f.readReal64();
+  h.time = f.readReal64();
+  h.redshift = f.readReal64();
+  h.flag_sfr = f.readInt32();
+  h.flag_feedback = f.readInt32();
+  for (int i = 0; i < 6; i++)
+    h.npartTotal[i] = f.readInt32();
+  h.flag_cooling = f.readInt32();
+  h.num_files = f.readInt32();
+  h.BoxSize = f.readReal64();
+  h.Omega0 = f.readReal64();
+  h.OmegaLambda = f.readReal64();
+  h.HubbleParam = f.readReal64();
+  f.endCheckpoint(true);
+  
+  data->NumPart = 0;
+  for(int k=0; k<5; k++)
+    data->NumPart += h.npart[k];
+  
+  data->pos = new FCoordinates[data->NumPart];
+  
+  f.beginCheckpoint();
+  for(int k = 0, p = 0; k < 5; k++) {
+    for(int n = 0; n < h.npart[k]; n++) {
+      data->pos[p][0] = f.readReal32();
+      data->pos[p][1] = f.readReal32();
+      data->pos[p][2] = f.readReal32();
+      p++;
+    }
+  }
+  f.endCheckpoint();
+  
+  // Skip velocities
+  f.skip((long)data->NumPart*3+2*4);
+  // Skip ids
+
+  return data;
+}

src/loadGadget.hpp

@@ -0,0 +1,12 @@
+#ifndef __COSMO_LOAD_GADGET_HPP
+#define __COSMO_LOAD_GADGET_HPP
+
+#include "load_data.hpp"
+
+namespace CosmoTool {
+  
+  PurePositionData *loadGadgetPosition(const char *fname);
+  
+};
+
+#endif