Initial import

AUTHOR

+Matthieu Schaller
+matthieu.schaller@durham.ac.uk

README

+Note that there is no guaranty that this module compiles, works, returns
+sensible results or allows you to do science. It is probably not bug free.
+Please always check the results you get make sense.
+
+Run:
+
+python setup.py install --user
+
+to create and install the package. 
+On the cosma machines you will need the following modules:
+
+gnu_comp/c4/4.8.1
+hdf5/1.8.9-mt

error.hpp

+#ifndef ERROR_HPP
+#define ERROR_HPP
+
+#include <iostream>
+#include <cstdlib>
+ 
+#define error(text)  do{std::cerr << "ERROR:" << __FILE__ << ':' << __func__ << "()" << __LINE__ << ": " << text << '\n';  std::abort();}while(0)
+
+#define debug(text) do{std::cerr << "DEBUG:" << __FILE__ << ':' << __func__ << "()" << __LINE__ << ": " << text << '\n';}while(0)
+
+#define debug_variable(text) do{std::cerr << "DEBUG:" << __FILE__ << ':' << __func__ << "()" << __LINE__ << ": " << #text << "=" << text << '\n';} while(0)
+
+#endif //ERROR_HPP

example.py

+import numpy
+import eagle
+import sys
+import os
+
+
+sim='/cosma5/data/Eagle/ScienceRuns/Planck1/L0100N1504/PE/Z0p10_W1p00_E_3p0_0p3_ALPHA1p0e6_rhogas1_reposlim3p0soft_100mgas_cosma/data/'
+tag='028_z000p000'
+
+length_EA = eagle.readArray("FOF", sim, tag, "FOF/GroupLength")
+#lengthType_EA = eagle.readArray("FOF", sim_EA, tag, "FOF/GroupLengthType")
+offset_EA = eagle.readArray("FOF", sim, tag, "FOF/GroupOffset")
+
+#print min(offset_EA), max(offset_EA)
+
+offset_2 = numpy.cumsum(length_EA)
+
+print offset_2
+
+offset_2 = numpy.insert(offset_2, 0, 0)
+
+print offset_2
+#print min(offset_EA), max(offset_EA)
+
+
+a[b] = c
+
+
+R_500 = eagle.readArray("SUBFIND_GROUP", sim, tag, "FOF/Group_R_Crit500", numThreads=64)
+R_2500 = eagle.readArray("SUBFIND_GROUP", sim, tag, "FOF/Group_R_Crit2500", numThreads=64)
+M_500 = eagle.readArray("SUBFIND_GROUP", sim, tag, "FOF/Group_M_Crit500", numThreads=64) * 1e10
+
+ratio = R_500 / R_2500
+
+massBins = numpy.logspace(8.1, 15.9, 40)
+ratio_count = numpy.zeros(numpy.size(massBins))
+ratio_mean = numpy.zeros(numpy.size(massBins))
+ratio2_mean = numpy.zeros(numpy.size(massBins))
+
+for i in range(numpy.size(R_500)):
+    if M_500[i] > 1e6 and R_500[i] > 0.001 and R_2500[i] > 0.001:
+        bin = -1
+        for j in range(numpy.size(massBins)-1):
+            if M_500[i] > massBins[j] and M_500[i] < massBins[j+1]:
+                bin = j
+                break
+
+        if bin != -1:
+            ratio_count[bin] = ratio_count[bin] + 1
+            ratio_mean[bin] = ratio_mean[bin] + ratio[i]
+            ratio2_mean[bin] = ratio2_mean[bin] + ratio[i]*ratio[i]
+
+
+
+for j in range(numpy.size(massBins)):
+    if ratio_count[j] != 0:
+        ratio_mean[j] = ratio_mean[j] / ratio_count[j]
+        ratio2_mean[j] = ratio2_mean[j] / ratio_count[j]
+    else:
+        ratio_mean[j] = 0
+        ratio2_mean[j] = 0
+
+
+ratio_sigma = numpy.sqrt(ratio2_mean - ratio_mean * ratio_mean)
+
+file = open("non_parametric_100_z0p0.dat",'w')
+file.write("# M_500 mean sigma count\n")
+for j in range(numpy.size(massBins)-1):
+    file.write("%e %f %f %i\n"%(10**(0.5*(numpy.log10(massBins[j])+numpy.log10(massBins[j+1]))), ratio_mean[j], ratio_sigma[j], ratio_count[j]))
+file.close()

examples/example.py

+import numpy
+import eagle
+import sys
+import os
+
+
+sim='/cosma5/data/Eagle/ScienceRuns/Planck1/L0100N1504/PE/Z0p10_W1p00_E_3p0_0p3_ALPHA1p0e6_rhogas1_reposlim3p0soft_100mgas_cosma/data/'
+tag='028_z000p000'
+
+length_EA = eagle.readArray("FOF", sim, tag, "FOF/GroupLength")
+#lengthType_EA = eagle.readArray("FOF", sim_EA, tag, "FOF/GroupLengthType")
+offset_EA = eagle.readArray("FOF", sim, tag, "FOF/GroupOffset")
+
+#print min(offset_EA), max(offset_EA)
+
+offset_2 = numpy.cumsum(length_EA)
+
+print offset_2
+
+offset_2 = numpy.insert(offset_2, 0, 0)
+
+print offset_2
+#print min(offset_EA), max(offset_EA)
+
+
+a[b] = c
+
+
+R_500 = eagle.readArray("SUBFIND_GROUP", sim, tag, "FOF/Group_R_Crit500", numThreads=64)
+R_2500 = eagle.readArray("SUBFIND_GROUP", sim, tag, "FOF/Group_R_Crit2500", numThreads=64)
+M_500 = eagle.readArray("SUBFIND_GROUP", sim, tag, "FOF/Group_M_Crit500", numThreads=64) * 1e10
+
+ratio = R_500 / R_2500
+
+massBins = numpy.logspace(8.1, 15.9, 40)
+ratio_count = numpy.zeros(numpy.size(massBins))
+ratio_mean = numpy.zeros(numpy.size(massBins))
+ratio2_mean = numpy.zeros(numpy.size(massBins))
+
+for i in range(numpy.size(R_500)):
+    if M_500[i] > 1e6 and R_500[i] > 0.001 and R_2500[i] > 0.001:
+        bin = -1
+        for j in range(numpy.size(massBins)-1):
+            if M_500[i] > massBins[j] and M_500[i] < massBins[j+1]:
+                bin = j
+                break
+
+        if bin != -1:
+            ratio_count[bin] = ratio_count[bin] + 1
+            ratio_mean[bin] = ratio_mean[bin] + ratio[i]
+            ratio2_mean[bin] = ratio2_mean[bin] + ratio[i]*ratio[i]
+
+
+
+for j in range(numpy.size(massBins)):
+    if ratio_count[j] != 0:
+        ratio_mean[j] = ratio_mean[j] / ratio_count[j]
+        ratio2_mean[j] = ratio2_mean[j] / ratio_count[j]
+    else:
+        ratio_mean[j] = 0
+        ratio2_mean[j] = 0
+
+
+ratio_sigma = numpy.sqrt(ratio2_mean - ratio_mean * ratio_mean)
+
+file = open("non_parametric_100_z0p0.dat",'w')
+file.write("# M_500 mean sigma count\n")
+for j in range(numpy.size(massBins)-1):
+    file.write("%e %f %f %i\n"%(10**(0.5*(numpy.log10(massBins[j])+numpy.log10(massBins[j+1]))), ratio_mean[j], ratio_sigma[j], ratio_count[j]))
+file.close()

examples/profile.py

+import numpy
+import eagle
+import sys
+import os
+
+
+sim='/cosma5/data/Eagle/ScienceRuns/Planck1/L0100N1504/PE/Z0p10_W1p00_E_3p0_0p3_ALPHA1p0e6_rhogas1_reposlim3p0soft_100mgas_cosma/data/'
+tag='028_z000p000'
+
+
+R_500 = eagle.readArray("SUBFIND_GROUP", sim, tag, "FOF/Group_R_Crit500", numThreads=64)
+R_2500 = eagle.readArray("SUBFIND_GROUP", sim, tag, "FOF/Group_R_Crit2500", numThreads=64)
+M_500 = eagle.readArray("SUBFIND_GROUP", sim, tag, "FOF/Group_M_Crit500", numThreads=64) * 1e10
+
+ratio = R_500 / R_2500
+
+massBins = numpy.logspace(8.1, 15.9, 40)
+ratio_count = numpy.zeros(numpy.size(massBins))
+ratio_mean = numpy.zeros(numpy.size(massBins))
+ratio2_mean = numpy.zeros(numpy.size(massBins))
+
+for i in range(numpy.size(R_500)):
+    if M_500[i] > 1e6 and R_500[i] > 0.001 and R_2500[i] > 0.001:
+        bin = -1
+        for j in range(numpy.size(massBins)-1):
+            if M_500[i] > massBins[j] and M_500[i] < massBins[j+1]:
+                bin = j
+                break
+
+        if bin != -1:
+            ratio_count[bin] = ratio_count[bin] + 1
+            ratio_mean[bin] = ratio_mean[bin] + ratio[i]
+            ratio2_mean[bin] = ratio2_mean[bin] + ratio[i]*ratio[i]
+
+
+
+for j in range(numpy.size(massBins)):
+    if ratio_count[j] != 0:
+        ratio_mean[j] = ratio_mean[j] / ratio_count[j]
+        ratio2_mean[j] = ratio2_mean[j] / ratio_count[j]
+    else:
+        ratio_mean[j] = 0
+        ratio2_mean[j] = 0
+
+
+ratio_sigma = numpy.sqrt(ratio2_mean - ratio_mean * ratio_mean)
+
+file = open("non_parametric_100_z0p0.dat",'w')
+file.write("# M_500 mean sigma count\n")
+for j in range(numpy.size(massBins)-1):
+    file.write("%e %f %f %i\n"%(10**(0.5*(numpy.log10(massBins[j])+numpy.log10(massBins[j+1]))), ratio_mean[j], ratio_sigma[j], ratio_count[j]))
+file.close()

profile.py

+import numpy
+import eagle
+import sys
+import os
+
+
+sim='/cosma5/data/Eagle/ScienceRuns/Planck1/L0100N1504/PE/Z0p10_W1p00_E_3p0_0p3_ALPHA1p0e6_rhogas1_reposlim3p0soft_100mgas_cosma/data/'
+tag='028_z000p000'
+
+
+R_500 = eagle.readArray("SUBFIND_GROUP", sim, tag, "FOF/Group_R_Crit500", numThreads=64)
+R_2500 = eagle.readArray("SUBFIND_GROUP", sim, tag, "FOF/Group_R_Crit2500", numThreads=64)
+M_500 = eagle.readArray("SUBFIND_GROUP", sim, tag, "FOF/Group_M_Crit500", numThreads=64) * 1e10
+
+ratio = R_500 / R_2500
+
+massBins = numpy.logspace(8.1, 15.9, 40)
+ratio_count = numpy.zeros(numpy.size(massBins))
+ratio_mean = numpy.zeros(numpy.size(massBins))
+ratio2_mean = numpy.zeros(numpy.size(massBins))
+
+for i in range(numpy.size(R_500)):
+    if M_500[i] > 1e6 and R_500[i] > 0.001 and R_2500[i] > 0.001:
+        bin = -1
+        for j in range(numpy.size(massBins)-1):
+            if M_500[i] > massBins[j] and M_500[i] < massBins[j+1]:
+                bin = j
+                break
+
+        if bin != -1:
+            ratio_count[bin] = ratio_count[bin] + 1
+            ratio_mean[bin] = ratio_mean[bin] + ratio[i]
+            ratio2_mean[bin] = ratio2_mean[bin] + ratio[i]*ratio[i]
+
+
+
+for j in range(numpy.size(massBins)):
+    if ratio_count[j] != 0:
+        ratio_mean[j] = ratio_mean[j] / ratio_count[j]
+        ratio2_mean[j] = ratio2_mean[j] / ratio_count[j]
+    else:
+        ratio_mean[j] = 0
+        ratio2_mean[j] = 0
+
+
+ratio_sigma = numpy.sqrt(ratio2_mean - ratio_mean * ratio_mean)
+
+file = open("non_parametric_100_z0p0.dat",'w')
+file.write("# M_500 mean sigma count\n")
+for j in range(numpy.size(massBins)-1):
+    file.write("%e %f %f %i\n"%(10**(0.5*(numpy.log10(massBins[j])+numpy.log10(massBins[j+1]))), ratio_mean[j], ratio_sigma[j], ratio_count[j]))
+file.close()

readArray.hpp

+#ifndef READ_ARRAY_HPP
+#define READ_ARRAY_HPP
+
+
+#include <vector>
+#include <string>
+#include <sstream>
+#include <iterator>
+#include <cmath>
+#include <cassert>
+#include <hdf5.h>
+#include <sys/time.h>
+#include <pthread.h>
+
+#include "vector.hpp"
+#include "tuple.hpp"
+#include "error.hpp"
+
+enum fileType{FOF,            //FOF group arrays
+	      FOF_PARTICLES,  //Particle info in FOF files
+	      SNIP_FOF,            //FOF group arrays
+	      SNIP_FOF_PARTICLES,  //Particle info in FOF files
+	      SNAP, 
+	      SNIP, 
+	      SUBFIND,        // Subfind sub-halo arrays
+	      SUBFIND_GROUP,  // Subfind halo arrays
+	      SUBFIND_IDS,
+	      SNIP_SUBFIND,        // Subfind sub-halo arrays
+	      SNIP_SUBFIND_GROUP,  // Subfind halo arrays
+	      SNIP_SUBFIND_IDS,
+	      PARTDATA,
+	      SNIP_PARTDATA
+};
+
+template<typename T> inline const hid_t hdf5Type();
+template<> inline const hid_t hdf5Type<char*>(){return H5T_C_S1; }
+template<> inline const hid_t hdf5Type<int>(){return H5T_NATIVE_INT; }
+template<> inline const hid_t hdf5Type<size_t>(){return H5T_NATIVE_UINT64; }
+template<> inline const hid_t hdf5Type<unsigned int>(){return H5T_NATIVE_UINT32; }
+template<> inline const hid_t hdf5Type<unsigned long long>(){return H5T_NATIVE_UINT64; }
+template<> inline const hid_t hdf5Type<float>(){return H5T_NATIVE_FLOAT; }
+template<> inline const hid_t hdf5Type<double>(){return H5T_NATIVE_DOUBLE; }
+template<> inline const hid_t hdf5Type<Tuple<6, int> >(){return H5T_NATIVE_INT; }
+template<> inline const hid_t hdf5Type<Tuple<6, unsigned int> >(){return H5T_NATIVE_UINT32; }
+template<> inline const hid_t hdf5Type<Tuple<6, unsigned long long> >(){return H5T_NATIVE_UINT64; }
+template<> inline const hid_t hdf5Type<Tuple<6, size_t> >(){return H5T_NATIVE_UINT64; }
+template<> inline const hid_t hdf5Type<Tuple<6, float> >(){return H5T_NATIVE_FLOAT; }
+template<> inline const hid_t hdf5Type<Tuple<6, double> >(){return H5T_NATIVE_DOUBLE; }
+
+/**
+ * \brief Returns the name of the file for a given directory, tag, fileNumber and type of file
+ */
+inline std::string fileName(fileType type, const std::string& dir, const std::string& tag, int i=0)
+{
+  std::ostringstream oss;
+  switch(type)
+    {
+    case FOF:
+    case FOF_PARTICLES:
+      oss << dir << "/groups_" << tag << "/group_tab_" << tag << "." << i << ".hdf5";
+      return oss.str();
+    case SNIP_FOF:
+    case SNIP_FOF_PARTICLES:
+      oss << dir << "/groups_snip_" << tag << "/group_snip_tab_" << tag << "." << i << ".hdf5";
+      return oss.str();
+    case SUBFIND:
+    case SUBFIND_GROUP:
+    case SUBFIND_IDS:
+      oss << dir << "/groups_" << tag << "/eagle_subfind_tab_" << tag << "." << i << ".hdf5";
+      return oss.str();
+    case SNIP_SUBFIND:
+    case SNIP_SUBFIND_GROUP:
+    case SNIP_SUBFIND_IDS:
+      oss << dir << "/groups_snip_" << tag << "/eagle_subfind_snip_tab_" << tag << "." << i << ".hdf5";
+      return oss.str();
+    case SNIP:
+      oss << dir << "/snipshot_" << tag << "/snip_" << tag << "." << i << ".hdf5";
+      return oss.str();
+    case SNAP:
+      oss << dir << "/snapshot_" << tag << "/snap_" << tag << "." << i << ".hdf5";
+      return oss.str();
+    case PARTDATA:
+      oss << dir << "/particledata_" << tag << "/eagle_subfind_particles_" << tag << "." << i << ".hdf5";
+      return oss.str();
+    case SNIP_PARTDATA:
+      oss << dir << "/particledata_snip_" << tag << "/eagle_subfind_snip_particles_" << tag << "." << i << ".hdf5";
+      return oss.str();
+    default:
+      error("Type of files not supported");
+    }
+}
+
+inline std::string numFilesFieldName(fileType type)
+{
+  switch(type)
+    {
+    case FOF:
+    case FOF_PARTICLES:
+    case SNIP_FOF:
+    case SNIP_FOF_PARTICLES:
+    case SUBFIND:
+    case SUBFIND_GROUP:
+    case SUBFIND_IDS:
+    case SNIP_SUBFIND:
+    case SNIP_SUBFIND_GROUP:
+    case SNIP_SUBFIND_IDS:
+      return "/Header/NTask";
+    case SNAP:
+    case SNIP:
+    case PARTDATA:
+    case SNIP_PARTDATA:
+      return "/Header/NumFilesPerSnapshot";
+    default:
+      error("Type of file unknown");
+    }
+}
+
+
+/**
+ * \brief Returns the header field containg the total number of elements
+ */
+inline std::string totElementFieldName(fileType type)
+{
+  switch(type)
+    {
+    case FOF_PARTICLES: 
+    case SNIP_FOF_PARTICLES: 
+      return "/Header/ToTNids";
+    case FOF: 
+    case SUBFIND_GROUP:
+    case SNIP_FOF: 
+    case SNIP_SUBFIND_GROUP:
+      return "/Header/TotNgroups";
+    case SUBFIND: 
+    case SNIP_SUBFIND: 
+      return "/Header/TotNsubgroups";
+    case SUBFIND_IDS:
+    case SNIP_SUBFIND_IDS:
+      return "/Header/TotNids";
+    case SNAP:
+    case SNIP:
+    case PARTDATA: 
+    case SNIP_PARTDATA: 
+      return "/Header/NumPart_Total";
+    default:
+      error("Type of file unknown");
+    }
+}
+
+/**
+ * \brief Returns the header field containg the number of elements in this file
+ */
+inline std::string elementFieldName(fileType type)
+{
+  switch(type)
+    {
+    case FOF_PARTICLES: 
+    case SUBFIND_IDS:
+    case SNIP_FOF_PARTICLES: 
+    case SNIP_SUBFIND_IDS:
+      return "/Header/Nids";
+    case FOF: 
+    case SUBFIND_GROUP:
+    case SNIP_FOF:
+    case SNIP_SUBFIND_GROUP:
+      return "/Header/Ngroups";
+    case SUBFIND: 
+    case SNIP_SUBFIND: 
+      return "/Header/Nsubgroups";
+    case SNAP:
+    case SNIP:
+    case PARTDATA: 
+    case SNIP_PARTDATA: 
+      return "/Header/NumPart_ThisFile";
+    default:
+      error("Type of file unknown");
+    }
+}
+
+/**
+ * \brief Returns the name of this fileType
+ */
+inline std::string fileTypeName(fileType type)
+{
+  switch(type)
+    {
+    case FOF: return "FOF";
+    case SNIP_FOF: return "FOF (snipshot)";
+    case FOF_PARTICLES: return "FOF (particle info)";
+    case SNIP_FOF_PARTICLES: return "FOF (particle info, snipshot)";
+    case SUBFIND: return "Subfind";
+    case SUBFIND_GROUP: return "Subfind (group tabs)";
+    case SUBFIND_IDS: return "Subfind (IDs)";
+    case SNIP_SUBFIND: return "Subfind (snipshot)";
+    case SNIP_SUBFIND_GROUP: return "Subfind (group tabs, snipshot)";
+    case SNIP_SUBFIND_IDS: return "Subfind (IDs, snipshot)";
+    case SNIP: return "Snipshot";
+    case SNAP: return "Snapshot";
+    case PARTDATA: return "ParticleData";
+    case SNIP_PARTDATA: return "ParticleData (snipshot)";
+    default:
+      error("Type of file unknown");
+    }
+}
+
+inline int getGadgetType(fileType type, const std::string& arrayName)
+{
+  switch(type)
+    {
+    case FOF:
+    case FOF_PARTICLES:
+    case SUBFIND:
+    case SUBFIND_GROUP:
+    case SUBFIND_IDS:
+    case SNIP_FOF:
+    case SNIP_FOF_PARTICLES:
+    case SNIP_SUBFIND:
+    case SNIP_SUBFIND_GROUP:
+    case SNIP_SUBFIND_IDS:
+      return 0;
+    case SNIP:
+    case SNAP:
+    case PARTDATA:
+    case SNIP_PARTDATA:
+      if(arrayName[0] != '/')
+	error("Ill-formed array name. Must start with a /");
+      return atoi(arrayName.substr(9,1).c_str());
+    default:
+      error("Type of file unknown");
+    }
+}
+
+
+/**
+ *\brief  Returns a vector of strings containing the different elements of a path of the form /aaa/bbb/ccc
+ */
+inline std::vector<std::string> tokenizePath(const std::string& path)
+{
+  std::vector<std::string> tokens;
+  const char delimiter='/';
+  std::string::size_type lastPos = path.find_first_not_of(delimiter, 0);
+  std::string::size_type pos     = path.find_first_of(delimiter, lastPos);
+  while (pos != std::string::npos || lastPos !=std::string::npos)
+    {
+      tokens.push_back(path.substr(lastPos, pos - lastPos));
+      lastPos = path.find_first_not_of(delimiter, pos);
+      pos = path.find_first_of(delimiter, lastPos);
+    }
+  return tokens;
+}
+
+/**
+ * \brief Returns the value of the attribute 'name' in an open HDF5 group.
+ */
+template<typename T>
+inline T readAttributeFromGroup(hid_t grp, const std::string& name)
+{
+  T data;
+
+  const hid_t h_attr = H5Aopen(grp, name.c_str(), H5P_DEFAULT);
+  if(h_attr < 0)
+    error("Impossible to open attribute");
+
+  const herr_t h_err = H5Aread(h_attr, hdf5Type<T>(), &data);
+  if(h_err < 0)
+    error("Impossible to read attribute");
+
+  H5Aclose(h_attr);
+  return data;
+}
+
+
+/**
+ * \brief Returns the value of the attribute 'name' in an open HDF5 file.
+ */
+template<typename T>
+T readAttribute(hid_t h_file, const std::string& name)
+{
+  const std::vector<std::string> tokens = tokenizePath(name);
+  std::vector<hid_t> groups(tokens.size(), 0);
+
+  groups[0] = h_file;
+  for(size_t i(1); i<groups.size(); ++i)
+    groups[i] = H5Gopen1(groups[i-1], (tokens[i-1]).c_str());
+    
+  const T attr = readAttributeFromGroup<T>(groups.back(), tokens.back().c_str());
+
+  for(size_t i(1); i<groups.size(); ++i)
+    H5Gclose(groups[i]);
+  return attr;
+}
+
+/**
+ * \brief Returns the value of the attribute 'name' in an HDF5 file given by its dir and tag.
+ */
+template<typename T>
+inline T readAttribute(fileType type, const std::string& dir, const std::string& tag, const std::string& name)
+{
+ 
+  const std::string fname=fileName(type, dir, tag);
+  const hid_t h_file = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+  const T attr = readAttribute<T>(h_file, name);
+  H5Fclose(h_file);
+  return attr;
+}
+
+
+
+/**
+ * \brief Returns the type of the attribute 'name' in an open HDF5 group.
+ */
+inline hid_t readAttributeTypeFromGroup(hid_t grp, const std::string& name)
+{
+  const hid_t h_attr = H5Aopen(grp, name.c_str(), H5P_DEFAULT);
+  if(h_attr < 0)
+    error("Impossible to open attribute");
+
+  const herr_t h_type = H5Aget_type(h_attr);
+  if(h_type < 0)
+    error("Impossible to read attribute type");
+
+  H5Aclose(h_attr);
+  return h_type;
+}
+
+
+/**
+ * \brief Returns the type of the attribute 'name' in an open HDF5 file.
+ */
+inline hid_t readAttributeType(hid_t h_file, const std::string& name)
+{
+  const std::vector<std::string> tokens = tokenizePath(name);
+  std::vector<hid_t> groups(tokens.size(), 0);
+
+  groups[0] = h_file;
+  for(size_t i(1); i<groups.size(); ++i)
+    groups[i] = H5Gopen1(groups[i-1], (tokens[i-1]).c_str());
+    
+  const hid_t attr = readAttributeTypeFromGroup(groups.back(), tokens.back().c_str());
+
+  for(size_t i(1); i<groups.size(); ++i)
+    H5Gclose(groups[i]);
+  return attr;
+}
+
+
+/**
+ * \brief Returns the HDF5 type of a given attribute 'name' in an HDF5 file given its dir and tag.
+ */
+inline hid_t readAttributeType(fileType type, const std::string& dir, const std::string& tag, const std::string& name)
+{
+  const std::string fname=fileName(type, dir, tag);
+  const hid_t h_file = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+  const hid_t h_type = readAttributeType(h_file, name);
+  H5Fclose(h_file);
+  return h_type;
+}
+
+/**
+ * \brief Returns the size of the attribute 'name' in an open HDF5 group.
+ */
+inline hid_t readAttributeSizeFromGroup(hid_t grp, const std::string& name)
+{
+  const hid_t h_attr = H5Aopen(grp, name.c_str(), H5P_DEFAULT);
+  if(h_attr < 0)
+    error("Impossible to open attribute");
+
+  const hid_t h_space = H5Aget_space(h_attr);
+  if(h_space < 0)
+    error("Impossible to open attribute space");
+
+  const htri_t isSimple = H5Sis_simple(h_space);
+  if(isSimple < 0)
+    error("Impossible to extract information about the attribute's dataspace");
+  else if(isSimple == 0)
+    error("Attribute's dataspace is not simple");
+    
+  hsize_t dims[1];
+  const int ndims = H5Sget_simple_extent_dims(h_space, dims, 0);
+  if(ndims < 0)
+    error("Attribute's dataspace is weird");
+  else if(ndims == 0)
+    dims[0] = 1;
+
+  H5Sclose(h_space);
+  H5Aclose(h_attr);
+  return dims[0];
+}
+
+
+/**
+ * \brief Returns the size of the attribute 'name' in an open HDF5 file.
+ */
+inline int readAttributeSize(hid_t h_file, const std::string& name)
+{
+  const std::vector<std::string> tokens = tokenizePath(name);
+  std::vector<hid_t> groups(tokens.size(), 0);
+
+  groups[0] = h_file;
+  for(size_t i(1); i<groups.size(); ++i)
+    groups[i] = H5Gopen1(groups[i-1], (tokens[i-1]).c_str());
+    
+  const int size = readAttributeSizeFromGroup(groups.back(), tokens.back().c_str());
+
+  for(size_t i(1); i<groups.size(); ++i)
+    H5Gclose(groups[i]);
+  return size;
+}
+
+/**
+ * \brief Returns the size of a given attribute 'name' in an HDF5 file given its dir and tag.
+ */
+inline int readAttributeSize(fileType type, const std::string& dir, const std::string& tag, const std::string& name)
+{
+  const std::string fname=fileName(type, dir, tag);
+  const hid_t h_file = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+  const int size = readAttributeSize(h_file, name);
+  H5Fclose(h_file);
+  return size;
+}
+
+
+/**
+ * \brief Returns the value of the attribute 'name' in an open HDF5 file.
+ */
+template<typename T>
+T readUnits(hid_t h_file, const std::string& name)
+{
+  const std::vector<std::string> tokens = tokenizePath(name);
+  std::vector<hid_t> groups(tokens.size(), 0);
+
+  groups[0] = h_file;
+  for(size_t i(1); i<groups.size()-1; ++i)
+    groups[i] = H5Gopen1(groups[i-1], (tokens[i-1]).c_str());
+    
+  const hid_t h_data = H5Dopen1(groups[groups.size()-2], tokens[groups.size()-2].c_str());
+  const T attr = readAttributeFromGroup<T>(h_data, tokens.back().c_str());
+  H5Dclose(h_data);
+
+  for(size_t i(1); i<groups.size()-1; ++i)
+    H5Gclose(groups[i]);
+  return attr;
+}
+
+
+
+/**
+ * \brief Returns the type of the attribute 'name' in an open HDF5 group.
+ */
+inline hid_t readArrayTypeFromGroup(hid_t grp, const std::string& name)
+{
+  const hid_t h_attr = H5Dopen(grp, name.c_str(), H5P_DEFAULT);
+  if(h_attr < 0)
+    error("Impossible to open array");
+
+  const herr_t h_type = H5Dget_type(h_attr);
+  if(h_type < 0)
+    error("Impossible to read array type");
+
+  H5Dclose(h_attr);
+  return h_type;
+}
+
+
+/**
+ * \brief Returns the type of the attribute 'name' in an open HDF5 file.
+ */
+inline hid_t readArrayType(hid_t h_file, const std::string& name)
+{
+  const std::vector<std::string> tokens = tokenizePath(name);
+  std::vector<hid_t> groups(tokens.size(), 0);
+
+  // for(int i(0); i<tokens.size();++i)
+  //   std::cout << tokens[i] << std::endl;
+
+  groups[0] = h_file;
+  for(size_t i(1); i<groups.size(); ++i)
+    {
+      // std::cout << "Opening group '" << tokens[i-1] << "'" <<std::endl;
+      groups[i] = H5Gopen1(groups[i-1], (tokens[i-1]).c_str());
+      // std::cout << "Group open" << std::endl;
+    }
+    
+  const hid_t attr = readArrayTypeFromGroup(groups.back(), tokens.back().c_str());
+
+  for(size_t i(1); i<groups.size(); ++i)
+    H5Gclose(groups[i]);
+  return attr;
+}
+
+
+/**
+ * \brief Returns the HDF5 type of a given attribute 'name' in an HDF5 file given its dir and tag.
+ */
+inline hid_t readArrayType(fileType type, const std::string& dir, const std::string& tag, const std::string& name)
+{
+  const std::string fname=fileName(type, dir, tag);
+  const hid_t h_file = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+  const hid_t h_type = readArrayType(h_file, name);
+  H5Fclose(h_file);
+  return h_type;
+}
+
+/**
+ * \brief Returns the size of the array 'name' in an open HDF5 group.
+ */
+inline hid_t readArraySizeFromGroup(hid_t grp, const std::string& name)
+{
+  const hid_t h_attr = H5Dopen(grp, name.c_str(), H5P_DEFAULT);
+  if(h_attr < 0)
+    error("Impossible to open dataset");
+
+  const hid_t h_space = H5Dget_space(h_attr);
+  if(h_space < 0)
+    error("Impossible to open data space");
+
+  const htri_t isSimple = H5Sis_simple(h_space);
+  if(isSimple < 0)
+    error("Impossible to extract information about the array's dataspace");
+  else if(isSimple == 0)
+    error("array's dataspace is not simple");
+    
+  hsize_t dims[2];
+  const int ndims = H5Sget_simple_extent_dims(h_space, dims, 0);
+  if(ndims < 0)
+    error("Array's dataspace is weird");
+
+  H5Sclose(h_space);
+  H5Dclose(h_attr);
+  if(ndims == 1)
+    return 1;
+  else if(ndims == 2)
+    return dims[1];
+  else
+    error("Array's dataspace is weird");
+}
+
+
+/**
+ * \brief Returns the size of the array 'name' in an open HDF5 file.
+ */
+inline int readArraySize(hid_t h_file, const std::string& name)
+{
+  const std::vector<std::string> tokens = tokenizePath(name);
+  std::vector<hid_t> groups(tokens.size(), 0);
+
+  groups[0] = h_file;
+  for(size_t i(1); i<groups.size(); ++i)
+    groups[i] = H5Gopen1(groups[i-1], (tokens[i-1]).c_str());
+    
+  const int size = readArraySizeFromGroup(groups.back(), tokens.back().c_str());
+
+  for(size_t i(1); i<groups.size(); ++i)
+    H5Gclose(groups[i]);
+  return size;
+}
+
+/**
+ * \brief Returns the size of a given array 'name' in an HDF5 file given its dir and tag.
+ */
+inline int readArraySize(fileType type, const std::string& dir, const std::string& tag, const std::string& name)
+{
+  const std::string fname=fileName(type, dir, tag);
+  const hid_t h_file = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+  const int size = readArraySize(h_file, name);
+  H5Fclose(h_file);
+  return size;
+}
+
+
+
+template<typename T, typename S> 
+void readArray(hid_t file, const std::string& name, S* data, size_t offset=0)
+{
+  const hid_t h_data = H5Dopen1(file, name.c_str());
+  if(h_data < 0)
+    error("Impossible to open file");
+
+  const hid_t h_space = H5Dget_space(h_data);
+  if(h_space < 0)
+    error("Impossible to open dataspace");
+
+  hsize_t dims[2], dummy[2];
+  int rank = H5Sget_simple_extent_dims(h_space, dims, dummy);
+
+  if(rank == 2)
+    offset *= dims[1];
+
+  const hid_t h_err = H5Dread(h_data, hdf5Type<T>(), H5S_ALL, H5S_ALL, H5P_DEFAULT, reinterpret_cast<T*>(data) + offset);
+  if(h_err < 0)
+    error("Impossible to read dataset");
+
+  H5Dclose(h_data);
+}
+
+
+/**
+ * \brief Reads the content of the dataset 'name' into the vector 'data' in an open HDF5 file.
+ */
+template<typename T, typename S> 
+inline void readArray(hid_t file, const std::string& name, std::vector<S>& data, size_t offset=0)
+{
+  readArray(file, name, &data[0], offset);
+}
+
+template <typename S>
+struct thread_info
+{
+  int id;
+  int numSnap;
+  int numThreads;
+  fileType type;
+  int partType;
+  std::string dir;
+  std::string tag;
+  std::string arrayName;
+  S* array;
+  size_t* sizes;
+  size_t* offsets;
+  bool oldSubfindFix;
+};
+
+template<typename S>
+void* thread_readSize(void *thread_data)
+{
+  thread_info<S> t_data = *reinterpret_cast<thread_info<S>*>( thread_data );
+
+  for(int i(t_data.id); i<t_data.numSnap; i+=t_data.numThreads)
+    {
+      //std::cerr << "Thread " << t_data.id << "/" << t_data.numThreads << " reading file " << i << "/" << t_data.numSnap  << std::endl;
+      const hid_t file = H5Fopen(fileName(t_data.type, t_data.dir, t_data.tag, i).c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+      if(file == 0)
+	error("Impossible to open file");
+
+      if(!t_data.oldSubfindFix)
+	t_data.sizes[i] = readAttribute<Tuple<6, int> >(file, elementFieldName(t_data.type))[t_data.partType];
+      else
+	{
+	  t_data.sizes[i] = readAttribute<Tuple<6, int> >(file, "/Header/NumPart_ThisFile")[2];
+	}
+
+      H5Fclose(file);
+    }
+  return 0;
+}
+
+template<typename T, typename S>
+void* thread_readData(void *thread_data)
+{
+  thread_info<S> t_data = *reinterpret_cast<thread_info<S>*>( thread_data );
+
+  for(int i(t_data.id); i<t_data.numSnap; i+=t_data.numThreads)
+    {
+      // std::cerr << "Thread " << t_data.id << "/" << t_data.numThreads << " reading file " << i << "/" << t_data.numSnap << " " << t_data.arrayName << " " <<  t_data.array << " " <<  t_data.offsets[i] << " " << t_data.type << std::endl;
+      const hid_t file = H5Fopen(fileName(t_data.type, t_data.dir, t_data.tag, i).c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+      if(file == 0)
+	error("Impossible to open file");
+      if(t_data.sizes[i] > 0)
+	{
+	  readArray<T>(file, t_data.arrayName, t_data.array, t_data.offsets[i]);
+	}
+      H5Fclose(file);
+
+      // std::cerr << "Done" << std::endl;
+    }
+
+  return 0;
+}
+
+/**
+ * \brief Reads the content of the dataset 'arrayName' into the vector 'array' for a given file (type, dir, tag)
+ */
+template<typename T, typename S>
+void readArray(fileType type, const std::string& dir, const std::string& tag, const std::string& arrayName, std::vector<S>& array, bool physicalUnits=true, bool noH=true, int numThreads=16, bool verbose=true)
+{
+  error("Strange function call");
+
+  /* Check that we have at least one thread */
+  if(numThreads < 1)
+    error("Invalid number of threads. Value must be >0.");
+
+ /* Open file #0 to read global information */  
+  const hid_t h_file = H5Fopen(fileName(type, dir, tag).c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+  if(h_file < 0)
+    error("Impossible to read file #0");
+
+  /* Decide which of the 6 values for the number of elements to read if we are dealing with a sn[ai]pshot */
+  const int partType = getGadgetType(type, arrayName);
+
+  /* Decide which of the 6 values for the number of elements to read if we are dealing with a sn[ai]pshot */
+  const int numSnap = readAttribute<int>(h_file, numFilesFieldName(type));
+  const size_t totNumElem = readAttribute<Tuple<6, int> >(h_file, totElementFieldName(type))[partType];
+  const double hubbleParam = readAttribute<double>(h_file, "/Header/HubbleParam");
+  const double a = readAttribute<double>(h_file, "/Header/ExpansionFactor");
+  H5Fclose(h_file);
+ 
+  if(verbose)
+    std::cout << "Extracting " << totNumElem << " values for '"<< arrayName << "' (PartType=" << partType << ") in " << numSnap << " " << fileTypeName(type) << " files with tag '" << tag << "' using " << numThreads << " thread(s)." << std::endl;
+  
+  /* Prepare array */
+  array.resize(totNumElem);
+  if(totNumElem == 0)
+    return;
+
+  /* Separate arrays for parallel reads */
+  std::vector<size_t> sizes(numSnap, 0);
+  std::vector<size_t> offsets(numSnap, 0);
+  
+  /* Prepare data for threads*/
+  pthread_t* threads = new pthread_t[numThreads];
+  thread_info<S>* thread_out = new thread_info<S>[numThreads];
+  for(int i(0); i<numThreads; ++i)
+    {
+      thread_out[i].id = i;
+      thread_out[i].numSnap = numSnap;
+      thread_out[i].numThreads = numThreads;
+      thread_out[i].type = type;
+      thread_out[i].partType = partType;
+      thread_out[i].dir = dir;
+      thread_out[i].tag = tag;
+      thread_out[i].arrayName = arrayName;
+      thread_out[i].array = &array[0];
+      thread_out[i].sizes = &sizes[0];
+      thread_out[i].offsets = &offsets[0];
+      thread_out[i].oldSubfindFix = false;
+    }
+
+  timeval begin, end;
+  gettimeofday(&begin, NULL);
+
+  if(verbose)
+    std::cout << "Reading array lengths" << std::flush;
+
+  /* Compute sizes for parallel reads using threads */
+  for(int i(0); i<numThreads; ++i)
+    {
+      const int ret = pthread_create(&threads[i], 0, thread_readSize<S>, &thread_out[i]);
+      if(ret != 0)
+	error("Error while launching thread");
+    }
+
+  for (int i(0); i < numThreads; ++i) 
+    {
+      const int ret = pthread_join(threads[i], 0);
+      if(ret != 0)
+	error("Error while stopping thread");
+    } 
+
+  gettimeofday(&end, NULL);  
+  if(verbose)
+    std::cout << " took " << ((end.tv_sec-begin.tv_sec)*1000000 + end.tv_usec-begin.tv_usec) / 1e6 << "s" << std::endl;
+
+
+  /* Do a scan to convert from size to offset */
+  for(int i(1); i<numSnap; ++i)
+    {
+      offsets[i] = offsets[i-1] + sizes[i-1];
+      //      std::cout << i << " " << offsets[i] << " " << sizes[i] << std::endl;
+    }
+
+
+  /* Read data from files using threads */
+  gettimeofday(&begin, NULL);
+  if(verbose)
+    std::cout << "Reading data" << std::flush;
+
+  for(int i(0); i<numThreads; ++i)
+    {
+      const int ret = pthread_create(&threads[i], 0, thread_readData<T,S>, &thread_out[i]);
+      if(ret != 0)
+	error("Error while launching thread");
+    }
+
+  for (int i(0); i < numThreads; ++i) 
+    {
+      const int ret = pthread_join(threads[i], 0);
+      if(ret != 0)
+	error("Error while stopping thread");
+    } 
+
+  gettimeofday(&end, NULL);  
+  if(verbose)
+    std::cout << " took " << ((end.tv_sec-begin.tv_sec)*1000000 + end.tv_usec-begin.tv_usec) / 1e6 << "s" << std::endl;
+
+  /* Find a non-empty file for unit conversion */
+  int nonEmpty = -1;
+  for(int i(0); i<numSnap; ++i)
+    {
+      if(sizes[i] > 0)
+	{
+	  nonEmpty = i;
+	  break;
+	}
+    }
+  //debug_variable(nonEmpty);
+
+  /* Read factors for unit conversion */
+  const hid_t h_file2 = H5Fopen(fileName(type, dir, tag, nonEmpty).c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+  const float a_exponant = readUnits<float>(h_file2, arrayName + "/aexp-scale-exponent");
+  const float h_exponant = readUnits<float>(h_file2, arrayName + "/h-scale-exponent");
+  H5Fclose(h_file2);
+
+  /* Apply comoving to physical units conversion */
+  if(physicalUnits)
+    {
+      if(a_exponant != 0.)
+	{
+	  if(verbose)
+	    std::cout << "Converting to physical units. (Multiplication by a^" << a_exponant << ", a=" << a  << ")" << std::endl;
+	  const double factor = std::pow(a, static_cast<double>(a_exponant));
+	  for(size_t i(0); i<totNumElem; ++i)
+	    array[i] *= factor;
+	}
+      else
+	{
+	  if(verbose)
+	    std::cout << "Converting to physical units. No conversion needed !" << std::endl;
+	}
+    }
+
+  /* Apply h units conversion */
+  if(noH)
+    {
+      if(h_exponant != 0.)
+	{
+	  if(verbose)
+	    std::cout << "Converting to h-free units. (Multiplication by h^" << h_exponant << ", h=" << hubbleParam << ")" << std::endl;
+	  const double factor = std::pow(hubbleParam, static_cast<double>(h_exponant));
+	  for(size_t i(0); i<totNumElem; ++i)
+	    array[i] *= factor;
+	}
+      else
+	{
+	  if(verbose)
+	    std::cout << "Converting to h-free units. No conversion needed !" << std::endl;
+	}
+    }
+
+  delete[] threads;
+  delete[] thread_out;
+}
+
+
+
+
+
+
+
+
+template<typename T, typename S>
+void readMultipleArrays(fileType type, const std::string& dir, const std::string& tag, const std::vector<std::string>& arrayName, std::vector<std::vector<S> >& array, bool physicalUnits=true, bool noH=true, int numThreads=16, bool verbose=true)
+{
+  /* Check that we have at least one thread */
+  if(numThreads < 1)
+    error("Invalid number of threads. Value must be >0.");
+
+  /* Open file #0 to read global information */  
+  const hid_t h_file = H5Fopen(fileName(type, dir, tag).c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+  if(h_file < 0)
+    error("Impossible to read file #0");
+
+  const size_t numArrays = arrayName.size();
+  if(numArrays != array.size())
+    error("'arrayName' and 'array' have different sizes");
+
+  /* Decide which of the 6 values for the number of elements to read if we are dealing with a sn[ai]pshot */
+  const int partType = getGadgetType(type, arrayName[0]);
+
+  /* Check that the partType is the same for all arrays */
+  for(size_t i(1); i<arrayName.size(); ++i)
+    if(partType != getGadgetType(type, arrayName[i]))
+      error("Not all arrays have the same particle type");
+
+  /* Read general information from file #0 */
+  const int numSnap = readAttribute<int>(h_file, numFilesFieldName(type));
+  const int totNumElem = readAttribute<Tuple<6, int> >(h_file, totElementFieldName(type))[partType];
+  const double hubbleParam = readAttribute<double>(h_file, "/Header/HubbleParam");
+  const double a = readAttribute<double>(h_file, "/Header/ExpansionFactor");
+  H5Fclose(h_file);
+
+  if(verbose)
+    std::cout << "Extracting " << totNumElem << " values for many arrays (PartType=" << partType << ") in " << numSnap << " " << fileTypeName(type) << " files with tag '" << tag << "' using " << numThreads << " thread(s)." << std::endl;
+  
+  /* Prepare array */
+  for(size_t i(0); i<numArrays; ++i)
+    array[i].resize(totNumElem);
+  if(totNumElem == 0)
+    return;
+
+  /* Separate arrays for parallel reads */
+  std::vector<size_t> sizes(numSnap, 0);
+  std::vector<size_t> offsets(numSnap, 0);
+  
+  /* Prepare data for threads*/
+  pthread_t* threads = new pthread_t[numThreads];
+  thread_info<S>* thread_out = new thread_info<S>[numThreads];
+  for(int i(0); i<numThreads; ++i)
+    {
+      thread_out[i].id = i;
+      thread_out[i].numSnap = numSnap;
+      thread_out[i].numThreads = numThreads;
+      thread_out[i].type = type;
+      thread_out[i].partType = partType;
+      thread_out[i].dir = dir;
+      thread_out[i].tag = tag;
+      thread_out[i].arrayName = arrayName[0];
+      thread_out[i].array = &array[0][0];
+      thread_out[i].sizes = &sizes[0];
+      thread_out[i].offsets = &offsets[0];
+      thread_out[i].oldSubfindFix = false;
+    }
+
+  timeval begin, end;
+  gettimeofday(&begin, NULL);
+
+  if(verbose)
+    std::cout << "Reading array lengths" << std::flush;
+
+  /* Compute sizes for parallel reads using threads */
+  for(int i(0); i<numThreads; ++i)
+    {
+      const int ret = pthread_create(&threads[i], 0, thread_readSize<S>, &thread_out[i]);
+      if(ret != 0)
+	error("Error while launching thread");
+    }
+
+  for (int i(0); i < numThreads; ++i) 
+    {
+      const int ret = pthread_join(threads[i], 0);
+      if(ret != 0)
+	error("Error while stopping thread");
+    } 
+
+  gettimeofday(&end, NULL);  
+  if(verbose)
+    std::cout << " took " << ((end.tv_sec-begin.tv_sec)*1000000 + end.tv_usec-begin.tv_usec) / 1e6 << "s" << std::endl;
+
+
+  /* Do a scan to convert from size to offset */
+  for(int i(1); i<numSnap; ++i)
+    offsets[i] = offsets[i-1] + sizes[i-1];
+
+  /* Loop over the arrays to read */
+  for(size_t j(0); j<numArrays; ++j)
+    {
+
+      for(int i(0); i<numThreads; ++i)
+	{
+	  thread_out[i].arrayName = arrayName[j];
+	  thread_out[i].array = &array[j][0];
+	}
+
+      /* Read data from files using threads */
+      gettimeofday(&begin, NULL);
+      if(verbose)
+	std::cout << "Reading data (" << arrayName[j] << ") " << std::flush;
+
+      for(int i(0); i<numThreads; ++i)
+	{
+	  const int ret = pthread_create(&threads[i], 0, thread_readData<T,S>, &thread_out[i]);
+	  if(ret != 0)
+	    error("Error while launching thread");
+	}
+
+      for (int i(0); i < numThreads; ++i) 
+	{
+	  const int ret = pthread_join(threads[i], 0);
+	  if(ret != 0)
+	    error("Error while stopping thread");
+	} 
+
+      gettimeofday(&end, NULL);  
+      if(verbose)
+	std::cout << " took " << ((end.tv_sec-begin.tv_sec)*1000000 + end.tv_usec-begin.tv_usec) / 1e6 << "s" << std::endl;
+      
+      /* Find a non-empty file for unit conversion */
+      int nonEmpty = -1;
+      for(int i(0); i<numSnap; ++i)
+	{
+	  if(sizes[i] > 0)
+	    {
+	      nonEmpty = i;
+	      break;
+	    }
+	}
+
+      /* Read factors for unit conversion */
+      const hid_t h_file2 = H5Fopen(fileName(type, dir, tag, nonEmpty).c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+      const float a_exponant = readUnits<float>(h_file2, arrayName[j] + "/aexp-scale-exponent");
+      const float h_exponant = readUnits<float>(h_file2, arrayName[j] + "/h-scale-exponent");
+      H5Fclose(h_file2);
+
+      /* Apply comoving to physical units conversion */
+      if(physicalUnits && a_exponant != 0.)
+	{
+	  const double factor = std::pow(a, static_cast<double>(a_exponant));
+	  for(int i(0); i<totNumElem; ++i)
+	    array[j][i] *= factor;
+	}
+      
+      /* Apply h units conversion */
+      if(noH && h_exponant != 0.)
+	{
+	  const double factor = std::pow(hubbleParam, static_cast<double>(h_exponant));
+	  for(int i(0); i<totNumElem; ++i)
+	    array[j][i] *= factor;
+	}
+    }
+
+  delete[] threads;
+  delete[] thread_out;
+}
+
+
+
+#endif //READ_ARRAY_HPP

readEagle.cpp

+#include <Python.h>
+#include <numpy/ndarrayobject.h>
+//#include <numpy/ndarraytypes.h>
+ 
+#include <string>
+#include <iostream>
+#include <cctype>
+
+#undef H5T_NATIVE_FLOAT_g
+
+#include "readArray.hpp"
+
+#pragma GCC diagnostic ignored "-Wwrite-strings"
+
+
+inline fileType getFileType(const std::string& type)
+{
+  if(type == "FOF")
+    return FOF;
+  else if(type == "FOF_PARTICLES")
+    return FOF_PARTICLES;
+  else if(type == "SNIP_FOF")
+    return SNIP_FOF;
+  else if(type == "SNIP_FOF_PARTICLES")
+    return SNIP_FOF_PARTICLES;
+  else if(type == "SNAP")
+    return SNAP;
+  else if(type == "SNAPSHOT")
+    return SNAP;
+  else if(type == "SNIP")
+    return SNIP;
+  else if(type == "SNIPSHOT")
+    return SNIP;
+  else if(type == "SUBFIND")
+    return SUBFIND;
+  else if(type == "SUBFIND_GROUP")
+    return SUBFIND_GROUP;
+  else if(type == "SUBFIND_PARTICLES")
+    return SUBFIND_IDS;
+  else if(type == "SUBFIND_IDS")
+    return SUBFIND_IDS;
+  else if(type == "SNIP_SUBFIND")
+    return SNIP_SUBFIND;
+  else if(type == "SNIP_SUBFIND_GROUP")
+    return SNIP_SUBFIND_GROUP;
+  else if(type == "SNIP_SUBFIND_PARTICLES")
+    return SNIP_SUBFIND_IDS;
+  else if(type == "SNIP_SUBFIND_IDS")
+    return SNIP_SUBFIND_IDS;
+  else if(type == "PARTDATA")
+    return PARTDATA;
+  else if(type == "SNIP_PARTDATA")
+    return SNIP_PARTDATA;
+  else
+    error("Invalid fileType provided.");
+} 
+
+template<typename T>inline const char* pyFormatter();
+template<>inline const char* pyFormatter<int>(){ return "i";}
+template<>inline const char* pyFormatter<unsigned int>(){ return "I";}
+template<>inline const char* pyFormatter<unsigned long long>(){ return "K";}
+template<>inline const char* pyFormatter<float>(){ return "f";}
+template<>inline const char* pyFormatter<double>(){ return "d";}
+
+template<typename T> inline const NPY_TYPES pyType();
+template<>inline const NPY_TYPES pyType<int>(){ return NPY_INT;}
+template<>inline const NPY_TYPES pyType<unsigned int>(){ return NPY_UINT;}
+template<>inline const NPY_TYPES pyType<unsigned long long>(){ return NPY_ULONGLONG;}
+template<>inline const NPY_TYPES pyType<float>(){ return NPY_FLOAT;}
+template<>inline const NPY_TYPES pyType<double>(){ return NPY_DOUBLE;}
+ 
+template<int size, typename T>
+inline PyObject* returnArrayAttribute(fileType ftype, const std::string& dir, const std::string& tag, const std::string& name)
+{
+  npy_intp dims[1] = {size};
+  PyArrayObject* array = (PyArrayObject *) PyArray_SimpleNew(1, dims,  pyType<T>());
+  const Tuple<size, T> ret = readAttribute<Tuple<size, T> >(ftype, dir, tag, name);
+  for(int i(0); i<size; ++i)
+    ((T*) (array->data))[i] = ret[i];
+  return (PyObject*) array;
+}
+
+template<typename T>
+PyObject* returnAttribute(fileType ftype, const std::string& dir, const std::string& tag, const std::string& name, int size)
+{
+  if(size == 1)
+    {
+      const T ret = readAttribute<T>(ftype, dir, tag, name);
+      return Py_BuildValue(pyFormatter<T>(), ret);
+    }
+  else
+    {
+      switch(size)
+      	{
+      	case 6: return returnArrayAttribute<6,T>(ftype, dir, tag, name);
+	default:
+	  error("Atrribute size not supported. Buying Matthieu a drink might help.");
+      	}
+    }
+  return 0;
+}
+
+
+PyObject* readAttributePython(PyObject* self, PyObject *args, PyObject *kwargs)
+{
+  char *c_type, *c_dir, *c_tag, *c_name;
+  static char* list[] = {"fileType", "directory", "tag", "attribute", 0};
+  if(!PyArg_ParseTupleAndKeywords(args, kwargs, "ssss", list, &c_type, &c_dir, &c_tag, &c_name ))
+    return 0;
+
+  /* Convert to C++ string types for convenience */
+  const std::string type=c_type;
+  const std::string dir=c_dir;
+  const std::string tag=c_tag;
+  const std::string name=c_name;
+
+  /* Convert the string to the internal enum value */
+  const fileType ftype = getFileType(type);
+
+  /* Recover the type of the attribute */
+  const hid_t htype = readAttributeType(ftype, dir, tag, name);
+
+  /* Recover the type of the attribute */
+  const int size = readAttributeSize(ftype, dir, tag, name);
+
+  /* Call the version of the function that corresponds to the type */
+  if(H5Tequal(htype, H5T_NATIVE_INT))
+    {
+      H5Tclose(htype);
+      return returnAttribute<int>(ftype, dir, tag, name, size);
+    }
+  else if(H5Tequal(htype, H5T_NATIVE_UINT32))
+    {
+      H5Tclose(htype);
+      return returnAttribute<unsigned int>(ftype, dir, tag, name, size);
+    }
+  else if(H5Tequal(htype, H5T_NATIVE_UINT64))
+    {
+      H5Tclose(htype);
+      return returnAttribute<unsigned long long>(ftype, dir, tag, name, size);
+    }
+  else if(H5Tequal(htype, H5T_NATIVE_FLOAT))
+    {
+      H5Tclose(htype);
+      return returnAttribute<float>(ftype, dir, tag, name, size);
+    }
+  else if(H5Tequal(htype, H5T_NATIVE_DOUBLE))
+    {
+      H5Tclose(htype);
+      return returnAttribute<double>(ftype, dir, tag, name, size);
+    }
+  else
+    {
+      error("Unsupported data type in HDF5 file. Buying Matthieu a drink might help.");
+    }
+
+  return 0;
+}
+
+
+
+// --------------------------------------------------------------------------------------------------------------
+
+
+template<typename T>
+PyObject* returnArray(fileType ftype, const std::string& dir, const std::string& tag, const std::string& name, bool physicalUnits, bool noH, bool useCGS, bool verbose, int numThreads, bool oldSubfindFix=false, bool owlsStyle=false)
+{
+ /* Open file #0 to read global information */  
+  const hid_t h_file = H5Fopen(fileName(ftype, dir, tag).c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+  if(h_file < 0)
+    error("Impossible to read file #0");
+
+  /* Decide which of the 6 values for the number of elements to read if we are dealing with a sn[ai]pshot */
+  const int partType = getGadgetType(ftype, name);
+
+  /* Decide which of the 6 values for the number of elements to read if we are dealing with a sn[ai]pshot */
+  const int numSnap = readAttribute<int>(h_file, numFilesFieldName(ftype));
+  const size_t totNumElem = readAttribute<Tuple<6, size_t> >(h_file, oldSubfindFix ? "/Header/NumPart_Total" : totElementFieldName(ftype))[oldSubfindFix ? 2 : partType];
+  const double hubbleParam = readAttribute<double>(h_file, "/Header/HubbleParam");
+  const double a = readAttribute<double>(h_file, "/Header/Time");
+  H5Fclose(h_file);
+
+  const int length = readArraySize(ftype, dir, tag, name);
+
+  if(verbose)
+    std::cout << "Extracting " << totNumElem << " values for '"<< name << "' (PartType=" << partType << ") in " << numSnap << " " << fileTypeName(ftype) << " files with tag '" << tag << "' using " << numThreads << " thread(s)" << ( oldSubfindFix ? " (with old subfind fix)." : "." ) << std::endl;
+
+  /* Create array */
+  PyArrayObject* array = 0;
+
+  if(length == 1)
+    {
+      npy_intp dims[1] = {static_cast<npy_intp>(totNumElem)};
+      array = (PyArrayObject *) PyArray_SimpleNew(1, dims,  pyType<T>());
+    }
+  else
+    {
+      npy_intp dims[2] = {static_cast<npy_intp>(totNumElem), length};
+      array = (PyArrayObject *) PyArray_SimpleNew(2, dims,  pyType<T>());
+    }
+
+
+  T* c_array = (T*) array->data;
+
+
+  /* Separate arrays for parallel reads */
+  std::vector<size_t> sizes(numSnap, 0);
+  std::vector<size_t> offsets(numSnap, 0);
+  
+  /* Prepare data for threads*/
+  pthread_t* threads = new pthread_t[numThreads];
+  thread_info<T>* thread_out = new thread_info<T>[numThreads];
+  for(int i(0); i<numThreads; ++i)
+    {
+      thread_out[i].id = i;
+      thread_out[i].numSnap = numSnap;
+      thread_out[i].numThreads = numThreads;
+      thread_out[i].type = ftype;
+      thread_out[i].partType = partType;
+      thread_out[i].dir = dir;
+      thread_out[i].tag = tag;
+      thread_out[i].arrayName = name;
+      thread_out[i].array = &c_array[0];
+      thread_out[i].sizes = &sizes[0];
+      thread_out[i].offsets = &offsets[0];
+      thread_out[i].oldSubfindFix = oldSubfindFix;
+    }
+  
+  timeval begin, end;
+  gettimeofday(&begin, NULL);
+
+  if(verbose)
+    std::cout << "Reading array lengths" << std::flush;
+  
+ 
+  /* Compute sizes for parallel reads using threads */
+  for(int i(0); i<numThreads; ++i)
+    {
+      const int ret = pthread_create(&threads[i], 0, thread_readSize<T>, &thread_out[i]);
+      if(ret != 0)
+	error("Error while launching thread");
+    }
+
+  for (int i(0); i < numThreads; ++i) 
+    {
+      const int ret = pthread_join(threads[i], 0);
+      if(ret != 0)
+	error("Error while stopping thread");
+    } 
+ 
+  gettimeofday(&end, NULL);  
+  if(verbose)
+    std::cout << " took " << ((end.tv_sec-begin.tv_sec)*1000000 + end.tv_usec-begin.tv_usec) / 1e6 << "s" << std::endl;
+
+  // for(int i(0); i<numSnap; ++i)
+  //   sizes[i] *= length;
+
+  /* Do a scan to convert from size to offset */
+  for(int i(1); i<numSnap; ++i)
+      offsets[i] = offsets[i-1] + sizes[i-1];
+
+  if( offsets[numSnap-1] + sizes[numSnap-1] != totNumElem )
+    error( "Sizes mismatch !" );
+ 
+  /* Read data from files using threads */
+  gettimeofday(&begin, NULL);
+  if(verbose) 
+    std::cout << "Reading data" << std::flush;
+
+  for(int i(0); i<numThreads; ++i)
+    {
+      int ret = -1;
+      switch(length)
+	{
+	case 1:
+	  ret = pthread_create(&threads[i], 0, thread_readData<T, T >, &thread_out[i]); break;
+	case 3:
+	  ret = pthread_create(&threads[i], 0, thread_readData<T, Vector<T> >, &thread_out[i]); break;
+	case 6:
+	  ret = pthread_create(&threads[i], 0, thread_readData<T, Tuple<6,T> >, &thread_out[i]); break;
+	default:
+	  error("Array size not supported. Talk to Matthieu...");
+	}
+      if(ret != 0)
+	error("Error while launching thread");
+    }
+
+  for (int i(0); i < numThreads; ++i) 
+    {
+      const int ret = pthread_join(threads[i], 0);
+      if(ret != 0)
+	error("Error while stopping thread");
+    } 
+
+  gettimeofday(&end, NULL);  
+  if(verbose)
+    std::cout << " took " << ((end.tv_sec-begin.tv_sec)*1000000 + end.tv_usec-begin.tv_usec) / 1e6 << "s" << std::endl;
+
+  /* Find a non-empty file for unit conversion */
+  int nonEmpty = -1;
+  for(int i(0); i<numSnap; ++i)
+    {
+      if(sizes[i] > 0)
+	{
+	  nonEmpty = i;
+	  break;
+	}
+    }
+
+  if(physicalUnits || noH || useCGS)
+    {
+
+      /* Read factors for unit conversion */
+      const hid_t h_file2 = H5Fopen(fileName(ftype, dir, tag, nonEmpty).c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+      const float a_exponant = readUnits<float>(h_file2, name + "/aexp-scale-exponent");
+      const float h_exponant = readUnits<float>(h_file2, name + "/h-scale-exponent");
+      const double cgsConversion = readUnits<double>(h_file2, name + "/CGSConversionFactor");
+      H5Fclose(h_file2);
+      
+      /* Apply comoving to physical units conversion */
+      if(physicalUnits)
+	{
+	  if(a_exponant != 0.)
+	    {
+	      if(verbose)
+		std::cout << "Converting to physical units. (Multiplication by a^" << a_exponant << ", a=" << a  << ")" << std::endl;
+	      const double factor = std::pow(a, static_cast<double>(a_exponant));
+	      for(size_t i(0); i<totNumElem*length; ++i)
+		c_array[i] *= factor;
+	    }
+	  else
+	    {
+	      if(verbose)
+		std::cout << "Converting to physical units. No conversion needed !" << std::endl;
+	    }
+	}
+      
+      /* Apply h units conversion */
+      if(noH)
+	{
+	  if(h_exponant != 0.)
+	    {
+	      if(verbose)
+		std::cout << "Converting to h-free units. (Multiplication by h^" << h_exponant << ", h=" << hubbleParam << ")" << std::endl;
+	      const double factor = std::pow(hubbleParam, static_cast<double>(h_exponant));
+	      for(size_t i(0); i<totNumElem*length; ++i)
+		c_array[i] *= factor;
+	    }
+	  else
+	    {
+	      if(verbose)
+		std::cout << "Converting to h-free units. No conversion needed !" << std::endl;
+	    }
+	}
+      
+      /* Apply CGS units conversion */
+      if(useCGS)
+	{
+	  if(cgsConversion != 0.)
+	    {
+	      if(verbose)
+		std::cout << "Converting to CGS units. (Multiplication by " << cgsConversion << ")" << std::endl;
+	      for(size_t i(0); i<totNumElem*length; ++i)
+		c_array[i] *= cgsConversion;
+	    }
+	  else
+	    {
+	      if(verbose)
+		std::cout << "Converting to CGS units. No conversion needed !" << std::endl;
+	    }
+	}
+    }
+
+
+  delete[] threads;
+  delete[] thread_out;
+
+
+  return (PyObject*)array;
+}
+
+
+PyObject* readArrayPython(PyObject* self, PyObject *args, PyObject *kwargs)
+{
+  char *c_type, *c_dir, *c_tag, *c_name;
+  bool physicalUnits = true , noH = true, useCGS = false, verbose = true, oldSubfindFix = false, owlsStyle = false;
+  int numThreads = 16;
+  static char* list[] = {"fileType", "directory", "tag", "attribute", "physicalUnits", "noH", "useCGS", "verbose", "numThreads", "oldSubfindFix", "OWLSstyle", 0};
+  if(!PyArg_ParseTupleAndKeywords(args, kwargs, "ssss|iiiiiii", list, &c_type, &c_dir, &c_tag, &c_name, &physicalUnits, &noH, &useCGS, &verbose, &numThreads, &oldSubfindFix, &owlsStyle ))
+    return 0;
+
+  /* Convert to C++ string types for convenience */
+  const std::string type=c_type;
+  const std::string dir=c_dir;
+  const std::string tag=c_tag;
+  const std::string name=c_name;
+
+  /* Check that we have at least one thread */
+  if(numThreads < 1)
+    error("Invalid number of threads. Value must be >0.");
+
+  /* Convert the string to the internal enum value */
+  const fileType ftype = getFileType(type);
+
+  /* Recover the type of the arrray */
+  const hid_t htype = readArrayType(ftype, dir, tag, name);
+
+  /* Check validity of fix */
+  if((oldSubfindFix && ftype != SUBFIND_IDS && ftype != SNIP_SUBFIND_IDS))
+    error("Can only use subfind fix when reading subfind IDs.");
+
+  /* Call the version of the function that corresponds to the type */
+  if(H5Tequal(htype, H5T_NATIVE_INT))
+    {
+      H5Tclose(htype);
+      return returnArray<int>(ftype, dir, tag, name, physicalUnits, noH, useCGS, verbose, numThreads, oldSubfindFix, owlsStyle);
+    }
+  else if(H5Tequal(htype, H5T_NATIVE_UINT32))
+    {
+      H5Tclose(htype);
+      return returnArray<unsigned int>(ftype, dir, tag, name, physicalUnits, noH, useCGS, verbose, numThreads, oldSubfindFix, owlsStyle);
+    }
+  else if(H5Tequal(htype, H5T_NATIVE_UINT64))
+    {
+      H5Tclose(htype);
+      return returnArray<unsigned long long>(ftype, dir, tag, name, physicalUnits, noH, useCGS, verbose, numThreads, oldSubfindFix, owlsStyle);
+    }
+  else if(H5Tequal(htype, H5T_NATIVE_FLOAT))
+    {
+      H5Tclose(htype);
+      return returnArray<float>(ftype, dir, tag, name, physicalUnits, noH, useCGS, verbose, numThreads, oldSubfindFix, owlsStyle);
+    }
+  else if(H5Tequal(htype, H5T_NATIVE_DOUBLE))
+    {
+      H5Tclose(htype);
+      return returnArray<double>(ftype, dir, tag, name, physicalUnits, noH, useCGS, verbose, numThreads, oldSubfindFix, owlsStyle);
+    }
+  else
+    {
+      error("Unsupported data type in HDF5 file. Buying Matthieu a drink might help.");
+    }
+
+
+  return 0;
+}
+
+// --------------------------------------------------------------------------------------------------------------
+
+
+// PyObject* example(PyObject* self, PyObject *args, PyObject *kwargs)
+// {
+//   int argument1, argument2;
+//   float argument3;
+//   static char* list[] = {"argument1", "argument2", "argument3", 0};
+//   if(!PyArg_ParseTupleAndKeywords(args, kwargs, "iif", list, &argument1, &argument2, &argument3))
+//     return 0;
+
+//   std::cout << "my first function !" << std::endl;
+
+//   npy_intp dims[2] = {3,3};
+//   PyArrayObject* array = (PyArrayObject *) PyArray_SimpleNew(2, dims,  NPY_DOUBLE);
+//   double* c_array = (double*) array->data;
+
+//   int i;
+//   for(i=0; i<9; ++i)
+//     c_array[i] = i+1;
+
+//   return (PyObject*) array;
+// }
+
+static PyMethodDef functions[]={
+  //  {"example", (PyCFunction)example, METH_VARARGS|METH_KEYWORDS, "One comment"},
+  {"readAttribute", (PyCFunction)readAttributePython, METH_VARARGS|METH_KEYWORDS, "Read one attribute from an EAGLE output file."},
+  {"readArray", (PyCFunction)readArrayPython, METH_VARARGS|METH_KEYWORDS, "Read one array from an EAGLE output file."},
+  {0, 0, 0, 0}
+};
+
+PyMODINIT_FUNC initeagle(void)
+{
+  Py_InitModule("eagle", functions);
+
+  import_array();
+}

setup.py

+#!/usr/bin/python
+
+from distutils.core import setup
+from distutils.core import Extension
+
+
+module=Extension('eagle', 
+                 include_dirs = ['/cosma/local/Python/2.7.3/lib/python2.7/site-packages/numpy/core/include'],
+                 libraries = ['hdf5'],
+                 sources = ['readEagle.cpp']
+                 )
+
+setup(name = 'Eagle',
+      version = '1.0',
+      description='Read EAGLE files in serial or parallel mode.',
+      ext_modules = [module],
+      author='Matthieu Schaller',
+      author_email='matthieu.schaller@durham.ac.uk',
+      url='http://eagle.strw.leidenuniv.nl/wiki/doku.php?id=eagle:documentation:reading_python',
+      long_description=''' See website. '''
+      )

tuple.hpp

+#ifndef TUPLE_HPP
+#define TUPLE_HPP
+
+#include <ostream>
+
+template<size_t N, typename T>
+struct Tuple
+{
+  const T& operator[](size_t i) const {return m_data[i];}
+  T& operator[](size_t i) {return m_data[i];}
+
+  template<typename S>
+  void operator*=(const S& other)
+  {
+    for(size_t i(0); i<N; ++i)
+      m_data[i] *= other;
+  }
+
+  size_t size() const
+  {
+    return N;
+  }
+
+  T m_data[N];
+
+};
+
+template<size_t N, typename T>
+std::ostream& operator<<(std::ostream& flux, const Tuple<N, T>& tuple)
+{
+  for(int i(0); i<N-1; ++i)
+    flux << tuple[i] << " ";
+  flux << tuple[N-1];
+  return flux;
+}
+
+template<size_t N, typename T>
+inline T sum(const Tuple<N, T>& tuple)
+{
+  T total(0);
+  for(size_t i(0); i<N; ++i)
+    total += tuple[i];
+  return total;
+}
+
+#endif //TUPLE_HPP

vector.hpp

+#ifndef VECTOR_HPP
+#define VECTOR_HPP
+
+/**
+ * \file Vector.hpp
+ * \author Matthieu Schaller
+ * \brief Defines a 3 dimensionnal vector class
+ * \date 16 December 2011
+ *
+ * Usual operators are overloaded and common mathematical functions are available
+ *
+ * \todo (Optimization) Implement meta-prog operators (Blitz++).
+ * \todo (Optimization) Implement SSE operations.
+ * \todo Write an optimized version of the inverseNorm function
+ */
+
+
+#include <cmath>
+#include <sstream>
+
+
+/**
+ * \class Vector<>
+ * \brief A 3D vector template class
+ *
+ * All the usual operators are overloaded. \n
+ * Components are in the public part of the class.
+ */
+template<typename T>
+struct Vector
+{
+
+  T x;  /**< X component*/
+  T y;  /**< Y component*/
+  T z;  /**< Z component*/
+
+  /**
+   *  \brief Builds a Vector with the 3 components specified
+   */
+  explicit Vector(const T& x_ = static_cast<T>(0),const T& y_ = static_cast<T>(0),const T& z_ = static_cast<T>(0))
+    :x(x_),y(y_),z(z_) {}
+
+
+  static const size_t dimension = 3; /**< Dimension of the vector*/
+
+  template<typename S>
+  Vector<T>& operator=(const Vector<S>& v)
+  {
+    x=v.x;
+    y=v.y;
+    z=v.z;
+    return *this;
+  }
+
+  /**
+   *  \brief Adds up two vectors
+   */
+  Vector<T>& operator+=(const Vector<T>& v)
+  {
+    x += v.x;
+    y += v.y;
+    z += v.z;
+    return *this;
+  }
+
+  /**
+   *  \brief Adds up a vector with a number
+   */
+  Vector<T>& operator+=(const T& alpha)
+  {
+    x += alpha;
+    y += alpha;
+    z += alpha;
+    return *this;
+  }
+
+  /**
+   *  \brief Substracts two vectors
+   */
+  Vector<T>& operator-=(const Vector<T>& v)
+  {
+    x -= v.x;
+    y -= v.y;
+    z -= v.z;
+    return *this;
+  }
+
+  /**
+   *  \brief Substracts a vector by a number
+   */
+  Vector<T>& operator-=(const T& alpha)
+  {
+    x -= alpha;
+    y -= alpha;
+    z -= alpha;
+    return *this;
+  }
+
+  /**
+   *  \brief Multiplies a vector by a number
+   */
+  Vector<T>& operator*=(const T& alpha)
+  {
+    x *= alpha;
+    y *= alpha;
+    z *= alpha;
+    return *this;
+  }
+  /**
+   *  \brief Divides a vector by a number
+   */
+  Vector<T>& operator/=(const T& alpha)
+  {
+    x /= alpha;
+    y /= alpha;
+    z /= alpha;
+    return *this;
+  }
+
+  /**
+   * \brief Increments all components by 1 and returns the new vector
+   */
+  Vector<T>& operator++ ()
+  {
+    ++x;
+    ++y;
+    ++z;
+    return *this;
+  }
+
+  /**
+   *  \brief Increments all components by 1 and returns the old vector
+   */
+  Vector<T> operator++ (int)
+  {
+    Vector<T> v(*this);
+    ++x;
+    ++y;
+    ++z;
+    return v;
+  }
+
+  /**
+   *  \brief Reduces all components by 1 and returns the new vector
+   */
+  Vector<T>& operator-- ()
+  {
+    --x;
+    --y;
+    --z;
+    return *this;
+  }
+
+  /**
+   *  \brief Reduces all components by 1 and return the old vector
+   */
+  Vector<T> operator-- (int)
+  {
+    Vector<T> v(*this);
+    --x;
+    --y;
+    --z;
+    return v;
+  }
+
+  /**
+   * \brief Access to the specified component
+   */
+  T& operator[](size_t i)
+  {
+    return *(reinterpret_cast<T*>(this) + i);         // Dirty trick. Don't do this at home
+  }
+
+  /**
+   * \brief Access to the specified component
+   */
+  const T& operator[](size_t i) const
+  {
+    return *(reinterpret_cast<const T*>(this) + i);         // Dirty trick. Don't do this at home
+  }
+
+};
+
+/**
+ *  \brief Compares two vectors component by component
+ */
+template<typename T>
+inline bool operator==(const Vector<T>& v,const Vector<T>& w)
+{
+  return (v.x == w.x && v.y == w.y && v.z == w.z);
+}
+
+/**
+ *  \brief Compares two vectors component by component
+ */
+template<typename T>
+inline bool operator!=(const Vector<T>& v,const Vector<T>& w)
+{
+  return (v.x != w.x || v.y != w.y || v.z != w.z);
+}
+
+/**
+ *  \brief Compares two vectors component by component
+ *  Vectors are sorted according to their x component then to their y component and finally z
+ *  This is useful for some operations implying STL containers and algorithms
+ */
+template<typename T>
+inline bool operator<(const Vector<T>& v,const Vector<T>& w)
+{
+  return ((v.x < w.x) ||
+	  (v.x == w.x && v.y < w.y) ||
+	  (v.x == w.x && v.y == w.y && v.z < w.z));
+}
+
+/**
+ *  \brief Compares two vectors component by component
+ */
+template<typename T>
+inline bool operator>=(const Vector<T>& v,const Vector<T>& w)
+{
+  return !(v < w);
+}
+
+/**
+ *  \brief Compares two vectors component by component
+ */
+template<typename T>
+inline bool operator>(const Vector<T>& v,const Vector<T>& w)
+{
+  return ((v.x > w.x) ||
+	  (v.x == w.x && v.y > w.y) ||
+	  (v.x == w.x && v.y == w.y && v.z > w.z));
+}
+
+/**
+ *  \brief Compares two vectors component by component
+ */
+template<typename T>
+inline bool operator<=(const Vector<T>& v,const Vector<T>& w)
+{
+  return !(v > w);
+}
+
+/**
+ *  \brief Returns the opposite of the Vector
+ */
+template<typename T>
+inline Vector<T> operator-(const Vector<T>& v)
+{
+  return Vector<T>(-v.x, -v.y, -v.z);
+}
+
+/**
+ *  \brief Applies the + operator to each of the Vector's component
+ */
+template<typename T>
+inline Vector<T> operator+(const Vector<T>& v)
+{
+  return Vector<T>(+v.x, +v.y, +v.z);
+}
+
+/**
+ *  \brief Computes the sum of two vectors
+ */
+template<typename T>
+inline Vector<T> operator+(const Vector<T>& v,const Vector<T>& w)
+{
+  return Vector<T>(v) += w;
+}
+
+/**
+ *  \brief Computes the difference between two vectors
+ */
+template<typename T>
+inline Vector<T> operator-(const Vector<T>& v,const Vector<T>& w)
+{
+  return Vector<T>(v) -= w;
+}
+
+/**
+ *  \brief Computes the product of a vector and a number
+ */
+template<typename T>
+inline Vector<T> operator*(const Vector<T>& v,const T& alpha)
+{
+  return Vector<T>(v) *= alpha;
+}
+
+/**
+ *  \brief Computes the product of a number and a vector
+ */
+template<typename T>
+inline Vector<T> operator*(const T& alpha,const Vector<T>& v)
+{
+  return Vector<T>(v) *= alpha;
+}
+
+/**
+ *  \brief Computes the ratio of a vector and a number
+ */
+template<typename T>
+inline Vector<T> operator/(const Vector<T>& v,const T& alpha)
+{
+  return Vector<T>(v) /= alpha;
+}
+
+/**
+ * \brief Computes the dot product of v and w
+ */
+template<typename T>
+inline T operator*(const Vector<T>& v,const Vector<T>& w)
+{
+  return v.x*w.x + v.y*w.y + v.z*w.z;
+}
+
+/**
+ * \brief Computes the dot product of v and w
+ */
+template<typename T>
+inline T dotProduct(const Vector<T>& v,const Vector<T>& w)
+{
+  return v.x*w.x + v.y*w.y + v.z*w.z;
+}
+
+/**
+ * \brief Computes the cross product of v and w
+ */
+template<typename T>
+inline Vector<T> crossProduct(const Vector<T>& v ,const Vector<T>& w)
+{
+  return Vector<T>(v.y * w.z - v.z * w.y,
+		   v.z * w.x - v.x * w.z,
+		   v.x * w.y - v.y * w.x);
+}
+
+/**
+ * \brief Computes the triple product of u,v and w
+ */
+template<typename T>
+T tripleProduct(const Vector<T>& u,const Vector<T>& v ,const Vector<T>& w)
+{
+  return u.x*v.y*w.z - u.x*v.z*w.y + u.y*v.z*w.x - u.y*v.x*w.z + u.z*v.x*w.y - u.z*v.y*w.x;
+}
+
+/**
+ * \brief Stream insertion operator
+ *
+ * Vectors are written as follow : (x,y,z)
+ */
+template<typename T, typename charT, class traits>
+std::basic_ostream<charT, traits>&
+operator<<(std::basic_ostream<charT, traits>& os, const Vector<T>& v)
+{
+  std::basic_ostringstream<charT, traits> s;
+  s.flags(os.flags());
+  s.imbue(os.getloc());
+  s.precision(os.precision());
+  s << v.x << " " << v.y << " " << v.z;
+  return os << s.str();
+}
+
+
+/**
+ * \brief Stream extraction operator
+ *
+ * Vectors of the follwing form will be accepted:\n
+ * (x,y,z) \n
+ * (x) \n
+ * x \n
+ * The latter two producing a vector whose components are (x,x,x)
+ * Any other input will raise the stream failbit.
+ *
+ */
+template<typename T, typename charT, class traits>
+std::basic_istream<charT, traits>&
+operator>>(std::basic_istream<charT, traits>& is, Vector<T>& v)
+{
+  T x, y, z;
+  charT temp;
+  is >> temp;
+  if (temp == '(')
+    {
+      is >> x >> temp;
+      if (temp == ',')
+        {
+	  is >> y >> temp;
+	  if (temp == ',')
+            {
+	      is >> z >> temp;
+	      if(temp == ')')
+		v = Vector<T>(x,y,z);
+	      else
+		is.setstate(std::ios_base::failbit);
+            }
+	  is.setstate(std::ios_base::failbit);
+        }
+      else if (temp == ')')
+	v = Vector<T>(x,x,x);
+      else
+	is.setstate(std::ios_base::failbit);
+    }
+  else
+    {
+      is.putback(temp);
+      is >> x;
+      v = Vector<T>(x,x,x);
+    }
+  return is;
+}
+
+/**
+ * \brief Computes the norm of the vector
+ */
+template <typename T>
+T norm(const Vector<T>& v)
+{
+  return std::sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
+}
+
+/**
+ * \brief Computes the square norm of the vector
+ */
+template <typename T>
+inline T norm2(const Vector<T>& v)
+{
+  return v.x*v.x + v.y*v.y + v.z*v.z;
+}
+
+/**
+ * \brief Computes the Manhattan norm of the vector
+ */
+template <typename T>
+inline T normManhattan(const Vector<T>& v)
+{
+  return std::abs(v.x) + std::abs(v.y) + std::abs(v.z);
+}
+
+/**
+ * \brief Computes inverse of the norm of the vector
+ */
+template<typename T>
+inline T inverseNorm(const Vector<T>& v)
+{
+  return 1. / norm(v);
+}
+
+/**
+ * \brief Computes cosinus of the angle between two vectors
+ */
+template <typename T>
+T cosAngle(const Vector<T>& v, const Vector<T>& w)
+{
+  return (v.x*w.x + v.y*w.y + v.z*w.z) / std::sqrt(norm2(v)*norm2(w));
+}
+
+/**
+ * \brief Computes the angle between two vectors
+ */
+template <typename T>
+T angle(const Vector<T>& v, const Vector<T>& w)
+{
+  return std::acos((v.x*w.x + v.y*w.y + v.z*w.z) / std::sqrt(norm2(v)*norm2(w)));
+}
+
+/**
+ * \brief Computes the projection of the vector b on the vector a
+ */
+template<typename T>
+Vector<T> proj(const Vector<T>& a,const Vector<T>& b)
+{
+  return Vector<T>(a)*=(a*b/norm2(a));
+}
+
+/**
+ * \brief Returns a unitary vector in the direction of a
+ */
+template<typename T>
+Vector<T> normalize(const Vector<T>& a)
+{
+  return Vector<T>(a)/norm(a);
+}
+
+/**
+ * \brief Returns a vector whose components are the absolute value of the vector
+ * given as an argument.
+ */
+template<typename T>
+Vector<T> abs(const Vector<T>& a)
+{
+  return Vector<T>(std::abs(a.x), std::abs(a.y), std::abs(a.z));
+}
+
+/**
+ * \brief Returns true if the vector is (0,0,0)
+ */
+template<typename T>
+inline bool isNull(const Vector<T>& v)
+{
+  return (v.x == static_cast<T>(0) && v.y == static_cast<T>(0) && v.z == static_cast<T>(0));
+}
+
+/**
+ * \brief Returns true if the vector is different from (0,0,0)
+ */
+template<typename T>
+inline bool isNotNull(const Vector<T>& v)
+{
+  return (v.x != static_cast<T>(0) || v.y != static_cast<T>(0) || v.z != static_cast<T>(0));
+}
+
+/**
+ * \brief Casts a vector of type T in a vector of type S. This is done by performing a static_cast on each component
+ */
+template<typename S,typename T>
+Vector<S> vector_cast(const Vector<T>& v)
+{
+  return Vector<S>(static_cast<S>(v.x),static_cast<S>(v.y),static_cast<S>(v.z));
+}
+
+
+
+#endif  // VECTOR_HPP