Imported Healpix, cfitsio, cosmotool. Added cmake tool to build dependencies (cfitsio, hdf5, netcdf, boost, healpix, gsl, ..). Adjusted CMakeLists.txt

external/cosmotool/src/algo.hpp

@@ -0,0 +1,63 @@
+#ifndef __COSMOTOOL_ALGO_HPP
+#define __COSMOTOOL_ALGO_HPP
+
+#include "config.hpp"
+
+namespace  CosmoTool
+{
+
+  template<typename T>
+  T cube(T a)
+  {
+    return a*a*a;
+  }
+
+  template<typename T>
+  T square(T a)
+  {
+    return a*a;
+  }
+
+  template<int N>
+  class SPowerBase
+  {
+  public:
+    template<typename T>
+    static T spower(T a)
+    {
+      if (N<0)
+        {
+          return 1/SPowerBase<-N>::spower(a);
+        }
+      if ((N%2)==0)
+	{
+	  T b = SPowerBase<N/2>::spower(a);
+	  return b*b;
+	}
+      T b = SPowerBase<(N-1)/2>::spower(a);
+      return a*b*b;
+    }
+  };
+
+  template<>
+  class SPowerBase<0>
+  {
+  public:
+    template<typename T>
+    static T spower(T a)
+    {
+      return T(1);
+    }
+  };
+
+  template<int N, typename T>
+  T spower(T a)
+  {
+    return SPowerBase<N>::spower(a);
+  }
+
+
+
+};
+
+#endif

external/cosmotool/src/bqueue.hpp

@@ -0,0 +1,30 @@
+#ifndef __COSMO_QUEUE_HPP
+#define __COSMO_QUEUE_HPP
+
+#include <cmath>
+
+namespace CosmoTool {
+
+  template<typename T,typename QType = int>
+  class BoundedQueue
+  {
+  public:
+    BoundedQueue(int maxSize, QType defaultMax);
+    BoundedQueue(T *pQueue, int maxSize, QType defaultMax);
+    ~BoundedQueue();
+
+    void push(T a, QType v);
+    T *getQueue() { return m_queue; }
+    QType *getPriorities() { return priority; }
+    QType getMaxPriority() { return priority[maxQueueSize-1]; }
+  private:
+    int maxQueueSize;
+    T *m_queue;
+    QType *priority;
+    bool autoFree;
+  };
+};
+
+#include "bqueue.tcc"
+
+#endif

external/cosmotool/src/bqueue.tcc

@@ -0,0 +1,67 @@
+namespace CosmoTool
+{
+  
+  template<typename T, typename QType>
+  BoundedQueue<T,QType>::BoundedQueue(int maxSize, QType defaultVal)
+  {
+    maxQueueSize = maxSize;
+    m_queue = new T[maxSize];
+    priority = new QType[maxSize];
+    autoFree = true;
+
+    for (int i = 0; i < maxSize; i++)
+      priority[i] = defaultVal;
+  }
+ 
+  template<typename T, typename QType>
+  BoundedQueue<T,QType>::BoundedQueue(T *pQueue, int maxSize, QType defaultVal)
+  {
+    maxQueueSize = maxSize;
+    m_queue = pQueue;
+    priority = new QType[maxSize];
+    autoFree = false;
+    
+    for (int i = 0; i < maxSize; i++)
+      priority[i] = defaultVal;
+
+    
+  }
+  
+  template<typename T, typename QType>
+  BoundedQueue<T,QType>::~BoundedQueue()
+  {
+    if (autoFree)
+      delete[] m_queue;
+    delete[] priority;
+  }
+
+
+  template<typename T, typename QType>
+  void BoundedQueue<T,QType>::push(T a, QType v)
+  {
+    if (v > priority[maxQueueSize-1])
+      return;
+
+    int i;
+    for (i = maxQueueSize-2; i >= 0; i--)
+      {
+	if (v > priority[i])
+	  {
+	    priority[i+1] = v;
+	    m_queue[i+1] = a;
+	    return;
+	  }
+	else
+	  {
+	    priority[i+1] = priority[i];
+	    m_queue[i+1] = m_queue[i];
+	  }
+      }
+    if (i < 0)
+      {
+	priority[0] = v;
+	m_queue[0] = a;
+      }
+  }
+   
+};

external/cosmotool/src/bsp_simple.hpp

@@ -0,0 +1,167 @@
+#ifndef __COSMOTOOL_SIMPLE_BSP_HPP
+#define __COSMOTOOL_SIMPLE_BSP_HPP
+
+#include <stdlib.h>
+#include <cmath>
+#include "algo.hpp"
+#include <queue>
+#include <exception>
+
+namespace CosmoTool
+{
+  
+  namespace simple_bsp
+  {
+
+    template<typename T, typename PType, int N>
+    struct space
+    {
+      typedef T data_t;
+      typedef PType point_t;
+      typedef PType coord_t[N];
+      static const int dim = N;
+    };
+
+    template<typename SType>
+    struct Plane
+    { 
+      typename SType::coord_t n;
+      typename SType::point_t d;
+    };
+
+    template<typename SType>
+    typename SType::point_t dot_product(const typename SType::coord_t& c1,
+				      const typename SType::coord_t& c2)
+    {
+      typename SType::point_t A = 0;
+
+      for(int j = 0; j < SType::dim; j++)
+	A += c1[j]*c2[j];
+      return A;
+    }
+
+    template<typename SType>
+    struct Node {
+      Plane<SType> plane;
+      Node<SType> *minus, *plus;
+      typename SType::data_t data;
+    };
+
+    template<typename SType>
+    void normal2(typename SType::coord_t p[2], typename SType::coord_t& n)
+    {
+      typename SType::point_t d;
+      using CosmoTool::square;
+
+      n[0] = p[1][1]-p[0][1];
+      n[1] = -p[1][0]+p[0][0];
+      d = std::sqrt(square(n[0])+square(n[1]));
+      n[0] /= d;
+      n[1] /= d;
+    }
+    
+    template<typename SType>
+    void normal3(typename SType::coord_t p[3], typename SType::coord_t& n)
+    {
+      typename SType::point_t delta0[3] = { p[1][0] - p[0][0], p[1][1] - p[0][1], p[1][2] - p[0][2] };
+      typename SType::point_t delta1[3] = { p[2][0] - p[0][0], p[2][1] - p[0][1], p[2][2] - p[0][2] };
+      typename SType::point_t d;
+      using CosmoTool::square;
+
+      n[0] = delta0[1] * delta1[2] - delta0[2] * delta1[1];
+      n[1] = delta0[2] * delta1[0] - delta0[0] * delta1[2];
+      n[2] = delta0[0] * delta1[1] - delta0[1] * delta1[0];
+      d = std::sqrt(square(n[0]) + square(n[1]) + square(n[2]));
+      n[0] /= d;
+      n[1] /= d;
+      n[2] /= d;
+    }
+    
+
+
+    template<typename SType>
+    struct Facet
+    {
+      typename SType::coord_t p[SType::dim];
+      typename SType::data_t data;
+
+      void center(typename SType::coord_t& c)
+      {
+	for (int j = 0; j < SType::dim; j++)
+	  {
+	    c[j] = 0;
+	    for (int i = 0; i < SType::dim; i++)
+	      {
+		c[j] += p[i][j]; 
+	      }
+	    c[j] /= SType::dim+1;
+	  }
+      }
+      
+      void normal(typename SType::coord_t& n)
+      {
+	if (SType::dim==2)
+	  {
+	    normal2<SType>(p, n);
+	    return;
+	  }
+	if (SType::dim == 3)
+	  {
+	    normal3<SType>(p, n);
+	    return;
+	  }
+	abort();
+      }
+
+    };
+
+
+    class InvalidPoint: public std::exception
+    {
+    };
+
+    template<typename T, typename PType, int N>
+    class BSP
+    {
+    public:
+      typedef space<T, PType, N> space_t;
+      typedef Plane<space_t> plane_t;
+      typedef Node<space_t> node_t;
+
+      node_t *root;
+      std::queue<node_t *> allocated;
+      
+      BSP() throw();
+      ~BSP();
+      void insert(Facet<space_t>& facet);
+
+      template<typename PIterator>
+      void insert(PIterator b, PIterator e, T data)
+      {
+	Facet<space_t> f;
+	int q = 0;
+
+	while (b != e && q < N+1)
+	  {
+	    for (int j = 0; j < N; j++)
+	      f.p[q][j] = (*b)[j];
+	    ++b;
+	    ++q;
+	  }
+	if (q != N)
+	  throw InvalidPoint();
+	
+	f.data = data;
+	insert(f);
+      }
+
+      bool inside(const typename space_t::coord_t& p) const;
+    };
+
+  };
+
+};
+
+#include "bsp_simple.tcc"
+
+#endif

external/cosmotool/src/bsp_simple.tcc

@@ -0,0 +1,117 @@
+#include <list>
+#include <queue>
+
+namespace CosmoTool
+{
+
+  namespace simple_bsp
+  {
+    
+    template<typename T, typename PType, int N>
+    BSP<T,PType,N>::BSP()
+     throw()
+    {
+      root = 0;
+    }
+
+    template<typename T, typename PType, int N>
+    BSP<T,PType,N>::~BSP()
+    {
+      while (!allocated.empty())
+	{
+	  node_t *r = allocated.front();
+	      
+	  allocated.pop();
+	  delete r;
+	}
+    }
+
+    template<typename T, typename PType, int N>
+    void BSP<T,PType,N>::insert(Facet<space_t>& f)
+    {
+      std::list<node_t **> subtrees;
+      std::list<node_t **> leaf_insert;
+
+      if (root != 0)
+	subtrees.push_back(&root);
+      else
+	leaf_insert.push_back(&root);
+
+      // Find the point of insertion. Do not bother to split triangle for this
+      // implementation.
+      while (!subtrees.empty())
+	{
+	  std::list<node_t **> new_subtrees;
+	  typename std::list<node_t **>::iterator iter = subtrees.begin();
+
+	  while (iter != subtrees.end())
+	    {
+	      typename space_t::point_t dp;
+	      bool cond_plus = false, cond_minus = false;
+	      node_t *current = *(*iter);
+	      
+	      for (int j = 0; j < N; j++)
+		{
+		  dp = dot_product<space_t>(f.p[j],current->plane.n) + current->plane.d;
+		  cond_plus = cond_plus || (dp > 0);
+		  cond_minus = cond_minus || (dp <= 0);
+		}
+	      
+	      bool joint = cond_plus && cond_minus;
+	      bool not_joint = (!cond_plus) && (!cond_minus);
+	      
+	      if (joint || not_joint)
+		{
+		  // Crawl and add another subtree
+		  *iter = &(current->minus);
+		  if (current->plus != 0)
+		    new_subtrees.push_back(&current->plus);
+		  else
+		    leaf_insert.push_back(&current->plus);
+		}
+	      else
+		{
+		  if (cond_plus)
+		    *iter = &current->plus;
+		  else
+		    *iter = &current->minus;
+		}
+	      if (*(*iter) == 0)
+		{
+		  leaf_insert.push_back(*iter);
+		  iter = subtrees.erase(iter);
+		}
+	      else
+		++iter;
+	    }
+	  if (!new_subtrees.empty())
+	    subtrees.splice(subtrees.end(), new_subtrees);
+	}
+      
+      node_t * current = new node_t;
+      f.normal(current->plane.n);
+      current->plane.d = -dot_product<space_t>((f.p[0]),current->plane.n);
+      
+      for (typename std::list<node_t **>::iterator i = leaf_insert.begin();
+	   i != leaf_insert.end();
+	   ++i)
+	*(*i) = current;
+
+      allocated.push(current);
+    }
+
+    template<typename T, typename PType, int N>
+    bool BSP<T,PType,N>::inside(const typename space_t::coord_t& p) const
+    {
+      node_t *current = root;
+
+      do
+	{
+	}
+      while();
+      current
+    }
+
+  };
+
+};

external/cosmotool/src/cic.cpp

@@ -0,0 +1,205 @@
+#include <assert.h>
+#include <math.h>
+#include <inttypes.h>
+#include "cic.hpp"
+
+CICFilter::CICFilter(uint32_t N, double len)
+{
+  spatialLen = len;
+  szGrid = N;
+  totalSize = N*N*N;
+  densityGrid = new CICType[totalSize];
+  resetMesh();
+}
+
+CICFilter::~CICFilter()
+{
+  delete[] densityGrid;
+}
+  
+void CICFilter::resetMesh()
+{
+  for (uint32_t i = 0; i < totalSize; i++)
+    densityGrid[i] = 0;
+}
+
+void CICFilter::putParticles(CICParticles *particles, uint32_t N)
+{
+#if 0
+  uint32_t numCorners = 1 << NUMDIMS;
+
+  for (uint32_t i = 0; i < N; i++)
+    {
+      Coordinates xyz;
+      int32_t ixyz[NUMDIMS];
+      int32_t rxyz[NUMDIMS];
+      CICType alpha[NUMDIMS];
+      CICType beta[NUMDIMS];
+      for (int j = 0; j < NUMDIMS; j++)
+	{
+	  xyz[j] = (particles[i].coords[j] / spatialLen * szGrid);
+	  ixyz[j] = (int32_t)floor(xyz[j] - 0.5);
+	  beta[j] = xyz[j] - ixyz[j] - 0.5;
+	  alpha[j] = 1 - beta[j];
+	  if (ixyz[j] < 0)
+	    ixyz[j] = szGrid-1;
+	}
+
+      CICType tot_mass = 0;
+      for (int j = 0; j < numCorners; j++)
+	{
+	  CICType rel_mass = 1;
+	  uint32_t idx = 0;
+	  uint32_t mul = 1;
+	  uint32_t mul2 = 1;
+
+	  for (int k = 0; k < NUMDIMS; k++)
+	    {
+	      uint32_t ipos = ((j & mul2) != 0);
+
+	      if (ipos == 1)
+		{
+		  rel_mass *= beta[k];
+		}
+	      else
+		{
+		  rel_mass *= alpha[k];
+		}
+
+	      rxyz[k] = ixyz[k] + ipos;
+	      
+	      if (rxyz[k] >= szGrid)
+		idx += (rxyz[k] - szGrid) * mul;
+	      else
+		idx += rxyz[k] * mul;
+
+	      mul2 *= 2;
+	      mul *= szGrid;
+	    }
+
+	  assert(rel_mass > 0);
+	  assert(rel_mass < 1);
+	  assert(idx < totalSize);
+	  densityGrid[idx] += rel_mass * particles[i].mass;
+	  tot_mass += rel_mass;
+	}
+      assert(tot_mass < 1.1);
+      assert(tot_mass > 0.9);
+    }
+#endif
+#if 0  
+  for (uint32_t i = 0; i < N; i++)
+    {
+      Coordinates xyz;
+      int32_t ixyz[NUMDIMS];
+      for (int j = 0; j < NUMDIMS; j++)
+	{
+	  xyz[j] = (particles[i].coords[j] / spatialLen * szGrid);
+	  ixyz[j] = (int32_t)round(xyz[j] - 0.5);
+	  if (ixyz[j] < 0)
+	    ixyz[j] = szGrid-1;
+	  else if (ixyz[j] >= szGrid)
+	    ixyz[j] = 0;
+	}
+
+      uint32_t idx = ixyz[0] + ixyz[1] * szGrid + ixyz[2] * szGrid * szGrid;
+      densityGrid[idx] += particles[i].mass;      
+    }
+
+#endif
+
+  for (uint32_t i = 0; i < N; i++)
+    {
+      CICType x, y, z;
+      int32_t ix, iy, iz;
+      int32_t ix2, iy2, iz2;
+
+      x = particles[i].coords[0] / spatialLen * szGrid + 0.5;
+      y = particles[i].coords[1] / spatialLen * szGrid + 0.5;
+      z = particles[i].coords[2] / spatialLen * szGrid + 0.5;
+
+      if (x < 0)
+	x += szGrid;
+      if (y < 0)
+	y += szGrid;
+      if (z < 0)
+	z += szGrid;
+
+      ix = ((int32_t)floor(x));
+      iy = ((int32_t)floor(y));
+      iz = ((int32_t)floor(z));
+      
+      ix2 = (ix + 1) % szGrid;
+      iy2 = (iy + 1) % szGrid;
+      iz2 = (iz + 1) % szGrid;
+
+      CICType alpha_x = x - ix;
+      CICType alpha_y = y - iy;
+      CICType alpha_z = z - iz;
+
+      ix %= szGrid;
+      iy %= szGrid;
+      iz %= szGrid;
+      
+      assert(alpha_x >= 0);
+      assert(alpha_y >= 0);
+      assert(alpha_z >= 0);
+
+      CICType beta_x = 1 - alpha_x;
+      CICType beta_y = 1 - alpha_y;
+      CICType beta_z = 1 - alpha_z;
+
+      assert(beta_x >= 0);
+      assert(beta_y >= 0);
+      assert(beta_z >= 0);
+
+      CICType mass = particles[i].mass;
+      uint32_t idx;
+
+      // 000
+      idx = ix + (iy + iz * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * beta_x * beta_y * beta_z;
+      
+      // 100
+      idx = ix2 + (iy + iz * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * alpha_x * beta_y * beta_z;
+      
+      // 010
+      idx = ix + (iy2 + iz * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * beta_x * alpha_y * beta_z;
+      
+      // 110
+      idx = ix2 + (iy2 + iz * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * alpha_x * alpha_y * beta_z;
+
+      // 001
+      idx = ix + (iy + iz2 * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * beta_x * beta_y * alpha_z;
+
+      // 101
+      idx = ix2 + (iy + iz2 * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * alpha_x * beta_y * alpha_z;
+
+      // 011
+      idx = ix + (iy2 + iz2 * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * beta_x * alpha_y * alpha_z;
+
+      // 111
+      idx = ix2 + (iy2 + iz2 * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * alpha_x * alpha_y * alpha_z;
+    }
+}
+  
+void CICFilter::getDensityField(CICType*& field, uint32_t& res)
+{
+  field = densityGrid;
+  res = totalSize;
+}

external/cosmotool/src/cic.hpp

@@ -0,0 +1,35 @@
+#ifndef __CICFILTER_HPP
+#define __CICFILTER_HPP
+
+#include "CosmoTool/config.hpp"
+#include <inttypes.h>
+
+using namespace CosmoTool;
+
+typedef float CICType;
+
+  typedef struct
+  {
+    float mass;
+    Coordinates coords;
+  } CICParticles;
+
+  class CICFilter
+  {
+  public:
+    CICFilter(uint32_t resolution, double spatialLen);
+    ~CICFilter();
+
+    void resetMesh();
+    void putParticles(CICParticles *particles, uint32_t N);    
+
+    void getDensityField(CICType*& field, uint32_t& res);
+
+  protected:
+    CICType *densityGrid;
+    double spatialLen;
+    uint32_t totalSize;
+    uint32_t szGrid;
+  };
+
+#endif

external/cosmotool/src/config.hpp

@@ -0,0 +1,134 @@
+#ifndef __COSMOTOOL_CONFIG_HPP
+#define __COSMOTOOL_CONFIG_HPP
+
+#include <string>
+#include <stdint.h>
+#include <exception>
+#include <cstring>
+
+namespace CosmoTool
+{
+
+#define NUMDIMS 3
+#define NUMCUBES 8
+
+  /**
+   * Base type to specity at what precision we
+   * must achieve computations.
+   */
+  typedef double ComputePrecision;
+  /**
+   * Coordinate type (should be a 3-array).
+   */
+  typedef double Coordinates[NUMDIMS];
+
+  /* 
+   * Single precision coordinates.
+   */
+  typedef float FCoordinates[NUMDIMS];
+  
+  /**
+   * This function is used whenever one needs a general
+   * conversion between mass and luminosity (or the opposite).
+   * It should take a "mass" (or luminosity) in input, a unit is
+   * given to convert this mass into solar units. The output should
+   * be the "luminosity" (or mass), in solar units.
+   */
+  typedef double (*BiasFunction)(double mass, double unit);
+  
+  /**
+   * Function to copy the coordinates "a" into "b".
+   */
+  inline void copyCoordinates(const Coordinates& a, Coordinates& b)
+  {
+    memcpy(b, a, sizeof(a));
+  }
+
+  /**
+   * Base exception class for all exceptions handled by
+   * this library.
+   */
+  class Exception : public std::exception
+  {
+  public:
+    Exception(const std::string& mess)
+      : msg(mess), msgok(true) {}
+    Exception()
+      : msgok(false) {}
+
+    virtual ~Exception() throw () {}
+
+    const char *getMessage() const { return msgok ? msg.c_str() : "No message"; };
+    virtual const char *what() const throw () { return msgok ? msg.c_str() : "What 'what' ?"; };
+
+  private:
+    std::string msg;
+    bool msgok;
+  };
+
+  /**
+   * Exception raised when an invalid argument has been
+   * passed to a function of the library.
+   */
+  class InvalidArgumentException : public Exception
+  {
+  public:
+    InvalidArgumentException(const std::string& mess)
+      : Exception(mess) {}
+    InvalidArgumentException()
+      : Exception() {}
+  };
+
+  /**
+   */
+  class InvalidRangeException : public Exception
+  {
+  public:
+    InvalidRangeException(const std::string& mess)
+      : Exception(mess) {}
+    InvalidRangeException()
+      : Exception() {}
+  };
+  
+  /**
+   */
+  class NoSuchFileException : public Exception
+  {
+  public:
+    NoSuchFileException(const std::string& mess)
+      : Exception(mess) {}
+    NoSuchFileException()
+      : Exception() {}
+  };
+
+  /**
+   */
+  class InvalidFileFormatException : public Exception
+  {
+  public:
+    InvalidFileFormatException(const std::string& mess)
+      : Exception(mess) {}
+    InvalidFileFormatException()
+      : Exception() {}
+  };
+
+  class EndOfFileException: public Exception
+  {
+  public:
+    EndOfFileException(const std::string& mess)
+      : Exception(mess) {}
+    EndOfFileException()
+      : Exception() {}
+  };
+
+  class FilesystemFullException: public Exception
+  {
+  public:
+    FilesystemFullException(const std::string& mess)
+      : Exception(mess) {}
+    FilesystemFullException()
+      : Exception() {}
+  };
+};
+
+#endif

external/cosmotool/src/dinterpolate.hpp

@@ -0,0 +1,94 @@
+#ifndef __COSMO_DINTERPOLATE_HPP
+#define __COSMO_DINTERPOLATE_HPP
+
+#include "config.hpp"
+#include "mykdtree.hpp"
+#include "kdtree_splitters.hpp"
+#include <gsl/gsl_eigen.h>
+
+namespace CosmoTool {
+ 
+  template<typename PType, typename IType, int N>
+  class DelaunayInterpolate {
+  public:
+    struct SimplexAccess {
+      int32_t *simplex_list;      
+    };
+ 
+    typedef KDTree<N, SimplexAccess, PType, KD_homogeneous_cell_splitter<N, SimplexAccess> > QuickTree;
+    typedef typename QuickTree::Cell QuickCell;
+    typedef PType CoordType[N];
+    
+    QuickTree *quickAccess;
+    QuickCell *cells;
+    PType *all_preweight;
+    int32_t *point_to_simplex_list_base;    
+    IType *values;
+    CoordType *positions;    
+    uint32_t numPoints;
+    uint32_t numSimplex;
+    uint32_t *simplex_list;
+    gsl_eigen_symmv_workspace *eigen_work;
+    bool *disable_simplex;
+    
+    /**
+     * This construct the interpolator. The construction is time consuming so
+     * please do it the less possible number of times, especially if you have
+     * a large number of points.
+     *
+     * @param positions list of the positions
+     * @param values list of the values taken at each position
+     * @param simplex_list list of points for each simplex. The packing
+     * is the following:
+     * [t(0,1),t(0,2),...,t(0,n+1),t(1,0),t(1,1),...,t(1,n+1),..],
+     * with t(i,j) the i-th simplex and j-th point of the simplex. The indexes
+     * refer to the previous list of points.
+     * @param numPoints the number of points
+     */
+    DelaunayInterpolate(CoordType *positions, IType *values, uint32_t *simplex_list,
+			uint32_t numPoints, uint32_t numSimplex)
+	throw (InvalidArgumentException)
+    {
+      this->positions = positions;
+      this->values = values;
+      this->simplex_list = simplex_list;
+      this->numPoints = numPoints;
+      this->numSimplex = numSimplex;
+      this->disable_simplex = new bool[numSimplex];
+      
+      buildPreweight();
+      buildQuickAccess();
+
+      eigen_work = gsl_eigen_symmv_alloc(N);
+
+    }
+    
+    ~DelaunayInterpolate()
+    {
+      delete[] cells;
+      delete quickAccess;
+      delete[] point_to_simplex_list_base;      
+      delete[] all_preweight;
+      delete[] disable_simplex;
+
+      gsl_eigen_symmv_free(eigen_work);
+    }
+
+    void buildPreweight()
+	throw (InvalidArgumentException);
+    void buildQuickAccess();
+    void buildHyperplane(const PType *v, CoordType& hyper);
+
+    bool checkPointInSimplex(const CoordType& pos, uint32_t simplex);
+    
+    uint32_t findSimplex(const CoordType& pos)
+      throw (InvalidArgumentException);
+
+    IType computeValue(const CoordType& pos)
+      throw (InvalidArgumentException);
+  };
+};
+
+#include "dinterpolate.tcc"
+
+#endif

external/cosmotool/src/dinterpolate.tcc

@@ -0,0 +1,339 @@
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <cstdlib>
+#include <cassert>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_linalg.h>
+
+namespace CosmoTool {
+
+
+  template<typename PType, typename IType, int N>
+  void DelaunayInterpolate<PType,IType,N>::buildQuickAccess()
+  {
+    cells = new QuickCell[numPoints];
+
+    uint32_t point_to_simplex_size = 0;
+    uint32_t *numSimplex_by_point = new uint32_t[numPoints];
+    uint32_t *index_by_point = new uint32_t[numPoints];
+
+    // First count the number of simplex for each point
+    for (uint32_t i = 0; i < numPoints; i++)
+      index_by_point[i] = numSimplex_by_point[i] = 0;
+    for (uint32_t i = 0; i < (N+1)*numSimplex; i++)
+      {
+	assert(simplex_list[i] < numPoints);
+	if (!disable_simplex[i/(N+1)])
+	  numSimplex_by_point[simplex_list[i]]++;
+      }
+
+    // Compute the total number and the index for accessing lists.
+    for (uint32_t i = 0; i < numPoints; i++)
+      {
+	index_by_point[i] = point_to_simplex_size;
+	point_to_simplex_size += numSimplex_by_point[i]+1;
+      }
+
+    // Now compute the real list.
+    point_to_simplex_list_base = new int32_t[point_to_simplex_size];
+    for (uint32_t i = 0; i < numSimplex; i++)
+      {
+	for (int j = 0; j <= N; j++)
+	  {
+	    uint32_t s = (N+1)*i+j;
+	    if (disable_simplex[i])
+	      continue;
+
+	    uint32_t p = simplex_list[s];
+	    assert(index_by_point[p] < point_to_simplex_size);
+	    point_to_simplex_list_base[index_by_point[p]] = i;
+	    ++index_by_point[p];
+	  }
+      }
+
+    // Finish the lists
+    for (uint32_t i = 0; i < numPoints; i++)
+      {
+	// check assertion
+	assert((i==0 && index_by_point[0]==numSimplex_by_point[0]) 
+	       || 
+	       ((index_by_point[i]-index_by_point[i-1]) == (numSimplex_by_point[i]+1)));
+	assert(index_by_point[i] < point_to_simplex_size);
+	point_to_simplex_list_base[index_by_point[i]] = -1;
+      }
+    
+    uint32_t idx = 0;
+    for (uint32_t i = 0; i < numPoints; i++)
+      {
+	cells[i].active = true;
+	cells[i].val.simplex_list = &point_to_simplex_list_base[idx];
+	// We may have to cast here.
+	for (int j = 0; j < N; j++)
+	  cells[i].coord[j] = positions[i][j];
+
+	idx += numSimplex_by_point[i]+1;
+      }
+
+    // Free the memory allocated for temporary arrays.
+    delete[] numSimplex_by_point;
+    delete[] index_by_point;
+
+    // Build the kd tree now.
+    quickAccess = new QuickTree(cells, numPoints);
+  }
+
+  template<typename PType, typename IType, int N>
+  void DelaunayInterpolate<PType,IType,N>::buildPreweight()
+	throw(InvalidArgumentException)
+  {
+    double preweight[N*N];
+    double preweight_inverse[N*N];
+    gsl_permutation *p = gsl_permutation_alloc(N);
+    uint32_t numDisabled = 0;
+
+    all_preweight = new PType[N*N*numSimplex];
+
+    for (uint32_t i = 0; i < numSimplex; i++)     
+      {
+	uint32_t base = i*(N+1);
+	uint32_t pref = simplex_list[base];
+	// Compute the forward matrix first.
+	for (int j = 0; j < N; j++)
+	  {
+	    PType xref = positions[pref][j];
+
+	    for (int k = 0; k < N; k++)
+	      {
+		preweight[j*N + k] = positions[simplex_list[k+base+1]][j] - xref;
+	      }
+	  }
+
+	gsl_matrix_view M = gsl_matrix_view_array(preweight, N, N);
+	gsl_matrix_view iM = gsl_matrix_view_array(preweight_inverse, N, N);
+	int signum;
+
+	gsl_linalg_LU_decomp(&M.matrix, p, &signum);
+        double a = fabs(gsl_linalg_LU_det(&M.matrix, signum));
+	if (a < 1e-10)
+         {
+#ifdef DEBUG
+           for (int j = 0; j < N; j++)
+            {
+              PType xref = positions[pref][j];
+
+              for (int k = 0; k < N; k++)
+                {
+                  preweight[j*N + k] = positions[simplex_list[k+base+1]][j] - xref;
+                }
+            }
+           std::ofstream f("matrix.txt");
+           for (int j = 0; j < N*N; j++)
+             f << std::setprecision(12) << preweight[j] << std::endl;
+	   throw InvalidArgumentException("Invalid tesselation. One simplex is coplanar.");
+#else
+           gsl_matrix_set_zero(&iM.matrix);
+	   disable_simplex[i] = true;
+	   numDisabled++;
+#endif
+         }
+        else {
+	  gsl_linalg_LU_invert(&M.matrix, p, &iM.matrix);
+	  disable_simplex[i] = false;
+        }
+ 
+	for (int j = 0; j < N*N; j++)
+	  all_preweight[N*N*i + j] = preweight_inverse[j];
+      }
+
+    std::cout << "Number of disabled simplices: " << numDisabled << std::endl;
+
+    gsl_permutation_free(p);
+  }
+
+  template<typename PType, typename IType, int N>
+  void DelaunayInterpolate<PType,IType,N>::buildHyperplane(const PType *v, CoordType& hyper)
+  {
+    double M[N][N], eVal[N], eVec[N][N];
+    gsl_matrix_view mM, evec;
+    gsl_vector_view eval;
+
+    // Construct the symmetric matrix
+    for (int k = 0; k < N; k++)
+      for (int l = k; l < N; l++)
+	{
+	  double val = 0;
+
+	  for (int i = 0; i < (N-1); i++)
+	    {
+	      val += v[i*N+l] * v[i*N+k];
+	    }
+	  M[l][k] = M[k][l] = val;
+	}
+
+    mM = gsl_matrix_view_array(&M[0][0], N, N);
+    evec = gsl_matrix_view_array(&eVec[0][0], N, N);
+    eval = gsl_vector_view_array(&eVal[0], N);
+
+    // Solve the eigensystem
+    gsl_eigen_symmv (&mM.matrix, &eval.vector, &evec.matrix, eigen_work);
+    
+    double minLambda = INFINITY;
+    uint32_t idx = N+1;
+
+    // Look for the smallest eigenvalue
+    for (int k = 0; k < N; k++)
+      {
+	if (minLambda > eVal[k])
+	  {
+	    minLambda = eVal[k];
+	    idx = k;
+	  }
+      }
+    assert(idx != (N+1));
+
+    // Copy the corresponding vector
+    for (int k = 0; k < N; k++)
+      {
+	hyper[k] = eVec[k][idx];
+      }
+  }
+
+  template<typename PType, typename IType, int N>
+  bool DelaunayInterpolate<PType,IType,N>::checkPointInSimplex(const CoordType& pos, uint32_t simplex)
+  {
+    if (disable_simplex[simplex])
+      return false;
+
+    uint32_t *desc_simplex = &simplex_list[simplex*(N+1)];
+    CoordType *p[N+1], v[N], hyper;
+
+    for (int k = 0; k <= N; k++)
+      p[k] = &positions[desc_simplex[k]];
+
+
+    for (int i = 0; i <= N; i++)
+      {
+	// Build vectors 
+	for (int k = 1; k <= N; k++)
+	  for (int l = 0; l < N; l++)
+	    v[k-1][l] = (*p[k])[l] - (*p[0])[l];
+
+	// Build hyperplane.
+	buildHyperplane(&v[0][0], hyper);
+       
+	// Compute the appropriate sign using the last point.
+	PType sign = 0;
+	for (int k = 0; k < N; k++)
+	  sign += hyper[k] * v[N-1][k];
+
+	// Now check the point has the same sign;
+	PType pnt_sign = 0;
+	for (int k = 0; k < N; k++)
+	  pnt_sign += hyper[k] * (pos[k] - (*p[0])[k]);
+
+	if (pnt_sign*sign < 0)
+	  return false;
+
+
+	// Rotate the points.
+	for (int k = 1; k <= N; k++)
+	  {
+	    p[k-1] = p[k];
+	  }
+	p[N] = &positions[desc_simplex[i]];
+      }
+
+    // We checked all possibilities. Return now.
+    return true;
+  }
+
+
+  template<typename PType, typename IType, int N>
+  uint32_t DelaunayInterpolate<PType,IType,N>::findSimplex(const CoordType& c)
+    throw (InvalidArgumentException)
+  {
+    uint32_t N_ngb = 1;
+    QuickCell **cell_Ngb = new QuickCell *[N_ngb];
+    typename QuickTree::coords kdc;
+
+    for (int i = 0; i < N; i++)
+       kdc[i] = c[i];
+
+    // It may happen that we are unlucky and have to iterate to farther
+    // neighbors. It is bound to happen, especially on the boundaries.
+    do
+      {
+	uint32_t i;
+
+	quickAccess->getNearestNeighbours(kdc, N_ngb, cell_Ngb);
+	
+	for (i = 0; i < N_ngb && cell_Ngb[i] != 0; i++)
+	  {
+	    int32_t *simplex_list = cell_Ngb[i]->val.simplex_list;
+	    uint32_t j = 0;
+
+	    while (simplex_list[j] >= 0)
+	      {
+		if (checkPointInSimplex(c, simplex_list[j]))
+		  {
+		    delete[] cell_Ngb;
+		    return simplex_list[j];
+		  }
+	        ++j;
+	      }
+	  }	
+	delete[] cell_Ngb;
+
+	// The point does not belong to any simplex.
+	if (i != N_ngb)
+	  throw InvalidArgumentException("the given point does not belong to any simplex");
+
+	N_ngb *= 2;
+	cell_Ngb = new QuickCell *[N_ngb];
+      }
+    while (1);
+
+    // Point not reached.
+    abort();
+    return 0;
+  }
+
+  template<typename PType, typename IType, int N>
+  IType DelaunayInterpolate<PType,IType,N>::computeValue(const CoordType& c)
+    throw (InvalidArgumentException)
+  {
+    uint32_t simplex = findSimplex(c);
+    PType *preweight = &all_preweight[simplex*N*N];
+    PType weight[N+1];
+    PType p0[N];
+    PType sum_weight = 0;
+
+    for (int i = 0; i < N; i++)
+      p0[i] = positions[simplex_list[simplex*(N+1) + 0]][i];
+
+    // Now we use the preweight to compute the weight...
+    for (int i = 1; i <= N; i++)
+      {
+	weight[i] = 0;
+	for (int j = 0; j < N; j++)
+	  weight[i] += preweight[(i-1)*N+j]*(c[j]-p0[j]);
+
+	assert(weight[i] > -1e-7);
+	assert(weight[i] < 1+1e-7);
+	sum_weight += weight[i];
+      }   
+    weight[0] = 1-sum_weight;
+    assert(weight[0] > -1e-7);
+    assert(weight[0] < (1+1e-7));
+    
+    // We compute the final value by weighing the value at the N+1
+    // points by the proper weight.
+    IType final = 0;
+    for (int i = 0; i <= N; i++)
+      final += weight[i] * values[ simplex_list[simplex*(N+1) + i] ];
+
+    return final;
+  }
+  
+};

external/cosmotool/src/eskow.hpp

@@ -0,0 +1,271 @@
+#ifndef __ESKOW_CHOLESKY_HPP
+#define __ESKOW_CHOLESKY_HPP
+
+#include <cmath>
+#include <vector>
+#include "mach.hpp"
+
+/* Implementation of Schnabel & Eskow, 1999, Vol. 9, No. 4, pp. 1135-148, SIAM J. OPTIM. */
+
+template<typename T, typename A>
+class CholeskyEskow
+{
+private:
+  static const bool verbose_eskow = true;
+  T tau, tau_bar, mu;
+
+  void print_matrix(A& m, int N)
+  {
+    using std::cout;
+    using std::endl;
+    using std::setprecision;
+
+    if (verbose_eskow)
+      {
+
+	for (int i = 0; i < N; i++)
+	  {
+	    for (int j = 0; j < N; j++)
+	      {
+		cout.width(6);
+		cout << setprecision(5) << m(i,j) << " ";
+	      }
+	    cout << endl;
+	  }
+	cout << endl;
+      }
+  }
+
+  T max_diag(A& m, int j, int N)
+  {
+    T maxval = std::abs(m(j,j));
+
+    for (int k = j+1; k < N; k++)
+      {
+	maxval = std::max(maxval, std::abs(m(k,k)));
+      }
+    return maxval;
+  }
+
+  void minmax_diag(A& m, int j, int N, T& minval, T& maxval, int& i_min, int& i_max)
+  {
+    i_min = i_max = j;
+    minval = maxval = m(j,j);
+
+    for (int k = j+1; k < N; k++)
+      {
+	maxval = std::max(maxval, m(k,k));
+	minval = std::min(minval, m(k,k));
+      }
+  
+    for (int k = j; k < N; k++)
+      {
+	if (m(k,k) == minval && i_min < 0)
+	  i_min = k;
+	if (m(k,k) == maxval && i_max < 0)
+	  i_max = k;
+      }
+  }
+
+  void swap_rows(A& m, int N, int i0, int i1)
+  {
+    for (int r = 0; r < N; r++)
+      std::swap(m(r,i0), m(r,i1));
+  }
+
+  void swap_cols(A& m, int N, int i0, int i1)
+  {
+    for (int c = 0; c < N; c++)
+      std::swap(m(i0,c), m(i1,c));
+  }
+
+  T square(T x)
+  {
+    return x*x;
+  }
+
+  T min_row(A& m, int j, int N)
+  {
+    T a = 1/m(j,j);
+    T v = m(j+1,j+1) - square(m(j+1,j))*a;
+    
+    for (int i = j+2; i < N; i++)
+      {
+	v = std::min(v, m(i, i) - square(m(i,j))*a);
+      }
+
+    return v;
+  }
+
+  int g_max(const std::vector<T>& g, int j, int N)
+  {
+    T a = g[j];
+    int k = j;
+
+    for (int i = j+1; i < N; i++)
+      {
+	if (a < g[i])
+	  {
+	    a = g[i];
+	    k = i;
+	  }
+      }
+    return k;
+  }
+
+public:
+  CholeskyEskow()
+  {
+    tau = std::pow(mach_epsilon<T>(), 1./3);
+    tau_bar = std::pow(mach_epsilon<T>(), 2./3);
+    mu=0.1;
+  }
+
+  void cholesky(A& m, int N, T& norm_E)
+  {
+    bool phaseone = true;
+    T gamma = max_diag(m, 0, N);
+    int j;
+
+    norm_E = 0;
+    
+    for (j = 0; j < N && phaseone; j++)
+      {
+	T minval, maxval;
+	int i_min, i_max;
+	
+	print_matrix(m, N);
+     
+	minmax_diag(m, j, N, minval, maxval, i_min, i_max);
+	if (maxval < tau_bar*gamma || minval < -mu*maxval)
+	  {
+	    phaseone = false;
+	    break;
+	  }
+      
+	if (i_max != j)	  
+	  {
+	    std::cout << "Have to swap i=" << i_max << " and j=" << j << std::endl;
+	    swap_cols(m, N, i_max, j);
+	    swap_rows(m, N, i_max, j);
+
+	  }
+      
+	if (min_row(m, j, N) < -mu*gamma)
+	  {
+	    phaseone = false;
+	    break;
+	  }
+      
+	T L_jj = std::sqrt(m(j,j));
+
+	m(j,j) = L_jj;
+	for (int i = j+1; i < N; i++)
+	  {
+	    m(i,j) /= L_jj;
+	    for (int k = j+1; k <= i; k++)
+	      m(i,k) -= m(i,j)*m(k,j);
+	  }
+      }
+
+
+    if (!phaseone && j == N-1)
+      {
+	T A_nn = m(N-1,N-1);
+	T delta = -A_nn + std::max(tau*(-A_nn)/(1-tau), tau_bar*gamma);
+      
+	m(N-1,N-1) = std::sqrt(m(N-1,N-1) + delta);	  	  
+      }
+
+    
+
+    if (!phaseone && j < (N-1))
+      {
+	std::cout << "Phase two ! (j=" << j << ")" << std::endl;
+
+	int k = j-1;
+	std::vector<T> g(N);
+
+	for (int i = k+1; i < N; i++)
+	  {
+	    g[i] = m(i,i);
+	    for (int j = k+1; j < i; j++)
+	      g[i] -= std::abs(m(i,j));
+	    for (int j = i+1; j < N; j++)
+	      g[i] -= std::abs(m(j,i));
+	  }
+
+	T delta, delta_prev = 0;
+	
+	for (int j = k+1; j < N-2; j++)
+	  {
+	    int i = g_max(g, j, N);
+	    T norm_j;
+
+	    print_matrix(m, N);
+
+	    if (i != j)
+	      {
+		swap_cols(m, N, i, j);
+		swap_rows(m, N, i, j);		
+	      }
+
+
+	    for (int i = j+1; j < N; j++)
+	      {
+		norm_j += std::abs(m(i,j));
+	      }
+	    
+	    delta = std::max(delta_prev, std::max((T)0, -m(j,j) + std::max(norm_j,tau_bar*gamma)));
+	    if (delta > 0)
+	      {
+		m(j,j) += delta;
+		delta_prev = delta;
+	      }
+	    
+	    if (m(j,j) != norm_j)
+	      {
+		T temp = 1 - norm_j/m(j,j);
+		
+		for (int i = j+1; j < N; j++)
+		  {
+		    g[i] += std::abs(m(i,j))*temp;
+		  }
+	      }
+
+	    // Now we do the classic cholesky iteration
+	    T L_jj = std::sqrt(m(j,j));
+	    
+	    m(j,j) = L_jj;
+	    for (int i = j+1; i < N; i++)
+	      {
+		m(i,j) /= L_jj;
+		for (int k = j+1; k <= i; k++)
+		  m(i,k) -= m(i,j)*m(k,j);
+	      }
+	  }
+
+	// The final 2x2 submatrix is special
+	T A00 = m(N-2, N-2), A01 = m(N-2, N-1), A11 = m(N-1,N-1);
+	T sq_DELTA = std::sqrt(square(A00-A11) + square(A01));
+	T lambda_hi = 0.5*((A00+A11) + sq_DELTA);
+	T lambda_lo = 0.5*((A00+A11) - sq_DELTA);
+       
+	delta = std::max(std::max((T)0, -lambda_lo + std::max(tau*sq_DELTA/(1-tau), tau_bar*gamma)),delta_prev);
+	if (delta > 0)
+	  {
+	    m(N-1,N-1) += delta;
+	    m(N,N) += delta;
+	    delta_prev = delta;
+	  }
+	m(N-2,N-2) = A00 = std::sqrt(A00);
+	m(N-1,N-2) = (A01 /= A00);
+	m(N-1,N-1) = std::sqrt(A11-A01*A01);
+	norm_E = delta_prev;
+      }
+  }
+
+};
+
+
+#endif

external/cosmotool/src/field.hpp

@@ -0,0 +1,83 @@
+#ifndef __COSMOTOOL_FIELD
+#define __COSMOTOOL_FIELD
+
+#include "config.hpp"
+#include <iostream>
+#include <cassert>
+
+namespace CosmoTool {
+
+  template<typename BaseType>
+  struct ScalarField
+  {
+    BaseType value;    
+  };
+
+  template<typename BaseType, int N>
+  struct VectorField
+  {
+    BaseType vec[N];
+    
+    VectorField& operator=(const VectorField& a)
+    {
+      for (int i = 0; i < N; i++)
+	vec[i] = a.vec[i];
+      return *this;
+    }
+
+    VectorField()
+    {
+      for (int i = 0; i < N; i++)
+	vec[i] = 0;      
+    }
+
+    VectorField(double a)
+    {
+      assert(a == 0);
+      for (int i = 0; i < N; i++)
+	vec[i] = 0;      
+    }
+  };
+
+  template<typename BaseType, int N>
+  VectorField<BaseType,N> operator*(BaseType s, const VectorField<BaseType,N>& a)
+  {
+    VectorField<BaseType,N> v;
+
+    for (int i = 0; i < N; i++)
+      v.vec[i] = a.vec[i]*s;
+    
+    return v;
+  }  
+  
+  template<typename BaseType, int N>
+  VectorField<BaseType,N> operator+(const VectorField<BaseType,N>& a, const VectorField<BaseType,N>& b)
+  {
+    VectorField<BaseType,N> v;
+
+    for (int i = 0; i < N; i++)
+      v.vec[i] = a.vec[i]+b.vec[i];
+    
+    return v;
+  }  
+
+  template<typename BaseType, int N>
+  VectorField<BaseType,N>& operator+=(VectorField<BaseType,N>& a, const VectorField<BaseType,N>& b)
+  {
+    for (int i = 0; i < N; i++)
+      a.vec[i]+=b.vec[i];
+    
+    return a;
+  }  
+
+};
+
+template<typename BaseType, int N>
+std::ostream& operator<<(std::ostream& s, const CosmoTool::VectorField<BaseType,N>& a)
+{
+  for (int i = 0; i < N; i++)
+    s << a.vec[i] << " " ;
+  return s;
+}
+
+#endif

external/cosmotool/src/fixArray.hpp

@@ -0,0 +1,40 @@
+#ifndef __FIX_ARRAY_HPP
+#define __FIX_ARRAY_HPP
+
+namespace CosmoTool
+{
+
+  template <typename T, unsigned int sz> class fixArray
+  {
+  private:
+    T d[sz];
+    
+  public:
+    /*! Returns the size of the array. */
+    long size() const { return sz; }
+
+    /*! Returns a reference to element \a #n */
+    template<typename T2> T &operator[] (T2 n) {
+      return d[n];
+    }
+
+    /*! Returns a constant reference to element \a #n */
+    template<typename T2> const T &operator[] (T2 n) const {
+      return d[n];
+    }
+
+    template<typename T2> void importArray(T2 *indata) {
+      for (int i = 0; i < sz; i++)
+	d[i] = indata[i];
+    }
+
+    template<typename T2> void exportArray(T2 *outdata) {
+      for (int i = 0; i < sz; i++)
+	outdata[i] = d[i];
+    }
+
+  };
+
+};
+
+#endif

external/cosmotool/src/fortran.cpp

@@ -0,0 +1,354 @@
+#include "config.hpp"
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include "fortran.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+UnformattedRead::UnformattedRead(const string& fname)
+  throw(NoSuchFileException)
+{
+  f = new ifstream(fname.c_str());
+  if (!*f)
+    throw NoSuchFileException();
+
+  swapOrdering = false;
+  cSize = Check_32bits;
+  checkPointRef = checkPointAccum = 0;
+}
+
+UnformattedRead::UnformattedRead(const char *fname)
+  throw(NoSuchFileException)
+{
+  f = new ifstream(fname);
+  if (!*f)
+    throw NoSuchFileException();
+  swapOrdering = false;
+  cSize = Check_32bits;
+  checkPointRef = checkPointAccum = 0;
+}
+
+
+UnformattedRead::~UnformattedRead()
+{
+  delete f;
+}
+
+// Todo implement primitive description
+void UnformattedRead::setOrdering(Ordering o)
+{
+  if (o == LittleEndian)
+    {
+    }
+}
+
+void UnformattedRead::setCheckpointSize(CheckpointSize cs)
+{
+  cSize = cs;
+}
+    
+void UnformattedRead::skip(int64_t off)
+  throw (InvalidUnformattedAccess)
+{
+  if (checkPointAccum == 0 && checkPointRef == 0)
+    {
+      // We are not in a checked block
+      f->seekg(off, ios::cur);
+      return;
+    }
+  if (off < 0)
+    throw InvalidUnformattedAccess();
+  
+  if ((checkPointAccum+off) > checkPointRef)
+       throw InvalidUnformattedAccess();
+
+  f->seekg(off, ios::cur);
+  checkPointAccum += off;
+}
+
+void UnformattedRead::beginCheckpoint()
+  throw (InvalidUnformattedAccess,EndOfFileException)
+{
+  if (checkPointAccum != 0)
+    throw InvalidUnformattedAccess();
+  
+  checkPointRef = (cSize == Check_32bits) ? 4 : 8;
+  checkPointAccum = 0;  
+
+  checkPointRef = (cSize == Check_32bits) ? readInt32() : readInt64();
+  checkPointAccum = 0;
+
+  if (f->eof())
+    throw EndOfFileException();
+}
+
+void UnformattedRead::endCheckpoint(bool autodrop)
+  throw (InvalidUnformattedAccess)
+{
+  if (checkPointRef != checkPointAccum)
+    {
+      if (!autodrop || checkPointAccum > checkPointRef)
+	throw InvalidUnformattedAccess();
+      f->seekg(checkPointRef-checkPointAccum, ios::cur);
+    }
+
+  int64_t oldCheckPoint = checkPointRef;
+
+  checkPointRef = (cSize == Check_32bits) ? 4 : 8;
+  checkPointAccum = 0;
+  checkPointRef = (cSize == Check_32bits) ? readInt32() : readInt64();
+  
+  if (oldCheckPoint != checkPointRef)
+    throw InvalidUnformattedAccess();
+
+  checkPointAccum = checkPointRef = 0;
+}
+
+void UnformattedRead::readOrderedBuffer(void *buffer, int size)
+      throw (InvalidUnformattedAccess)
+{
+  if ((checkPointAccum+size) > checkPointRef)
+       throw InvalidUnformattedAccess();
+
+  f->read((char *)buffer, size);
+  
+  if (swapOrdering)
+    {
+      char *cb = (char *)buffer;
+      for (int i = 0; i < size/2; i++)
+	swap(cb[i], cb[size-i-1]);
+    }
+  checkPointAccum += size;
+}
+  
+double UnformattedRead::readReal64()
+  throw (InvalidUnformattedAccess)
+{
+  union
+  {
+    char b[8];
+    double d;
+  } a;
+
+  readOrderedBuffer(&a, 8);
+
+  return a.d;
+}
+
+float UnformattedRead::readReal32()
+  throw (InvalidUnformattedAccess)
+{
+  union
+  {
+    char b[4];
+    float f;
+  } a;
+
+  readOrderedBuffer(&a, 4);
+
+  return a.f;
+}
+
+int32_t UnformattedRead::readInt32()
+  throw (InvalidUnformattedAccess)
+{
+  union
+  {
+    char b[4];
+    int32_t i;
+  } a;
+
+  readOrderedBuffer(&a, 4);
+
+  return a.i;
+}
+
+int64_t UnformattedRead::readInt64()
+  throw (InvalidUnformattedAccess)
+{
+  union
+  {
+    char b[8];
+    int64_t i;
+  } a;
+
+  readOrderedBuffer(&a, 8);
+
+  return a.i;
+}
+
+//// UnformattedWrite
+
+UnformattedWrite::UnformattedWrite(const string& fname)
+  throw(NoSuchFileException)
+{
+  f = new ofstream(fname.c_str());
+  if (!*f)
+    throw NoSuchFileException();
+
+  swapOrdering = false;
+  cSize = Check_32bits;
+  checkPointRef = checkPointAccum = 0;
+}
+
+UnformattedWrite::UnformattedWrite(const char *fname)
+  throw(NoSuchFileException)
+{
+  f = new ofstream(fname);
+  if (!*f)
+    throw NoSuchFileException();
+  swapOrdering = false;
+  cSize = Check_32bits;
+  checkPointRef = checkPointAccum = 0;
+}
+
+
+UnformattedWrite::~UnformattedWrite()
+{
+  delete f;
+}
+
+// Todo implement primitive description
+void UnformattedWrite::setOrdering(Ordering o)
+{
+  if (o == LittleEndian)
+    {
+    }
+}
+
+void UnformattedWrite::setCheckpointSize(CheckpointSize cs)
+{
+  cSize = cs;
+}
+    
+void UnformattedWrite::beginCheckpoint()
+  throw (InvalidUnformattedAccess,FilesystemFullException)
+{
+  if (checkPointAccum != 0)
+    throw InvalidUnformattedAccess();
+  
+  checkPointRef = f->tellp();
+  if (cSize == Check_32bits)
+    writeInt32(0);
+  else
+    writeInt64(0);
+
+  checkPointAccum = 0;
+
+  if (!*f)
+    throw FilesystemFullException();
+}
+
+void UnformattedWrite::endCheckpoint()
+  throw (InvalidUnformattedAccess,FilesystemFullException)
+{
+  if (checkPointAccum == 0)
+    throw InvalidUnformattedAccess();
+
+  streampos curPos = f->tellp();
+
+  int64_t deltaPos = curPos-checkPointRef;
+
+  deltaPos -= (cSize == Check_32bits) ? 4 : 8;
+  // This is a sanity check.
+  if (checkPointAccum != deltaPos)
+    throw InvalidUnformattedAccess();
+
+  uint64_t saveAccum = checkPointAccum;
+
+  f->seekp(checkPointRef);
+  if (cSize == Check_32bits)
+    writeInt32(saveAccum);
+  else
+    writeInt64(saveAccum);
+
+  f->seekp(curPos);
+  if (cSize == Check_32bits)
+    writeInt32(saveAccum);
+  else
+    writeInt64(saveAccum);
+
+  checkPointAccum = checkPointRef = 0;
+}
+
+void UnformattedWrite::writeOrderedBuffer(void *buffer, int size)
+  throw (FilesystemFullException)
+{
+  f->write((char *)buffer, size);
+  
+  if (swapOrdering)
+    {
+      char *cb = (char *)buffer;
+      for (int i = 0; i < size/2; i++)
+	swap(cb[i], cb[size-i-1]);
+    }
+  checkPointAccum += size;
+
+  if (!*f)
+    throw FilesystemFullException();
+}
+  
+void UnformattedWrite::writeReal64(double d)
+  throw (FilesystemFullException)
+{
+  union
+  {
+    char b[8];
+    double d;
+  } a;
+
+  a.d = d;
+
+  writeOrderedBuffer(&a, 8);
+}
+
+void UnformattedWrite::writeReal32(float f)
+  throw (FilesystemFullException)
+{
+  union
+  {
+    char b[4];
+    float f;
+  } a;
+  
+  a.f = f;
+
+  writeOrderedBuffer(&a, 4);
+}
+
+void UnformattedWrite::writeInt32(int32_t i)
+  throw (FilesystemFullException)
+{
+  union
+  {
+    char b[4];
+    int32_t i;
+  } a;
+
+  a.i = i;
+  writeOrderedBuffer(&a, 4);
+}
+
+void UnformattedWrite::writeInt64(int64_t i)
+  throw (FilesystemFullException)
+{
+  union
+  {
+    char b[8];
+    int64_t i;
+  } a;
+
+  a.i = i;
+  writeOrderedBuffer(&a, 8);
+}
+
+void UnformattedWrite::writeInt8(int8_t i)
+  throw (FilesystemFullException)
+{
+  union { char b; int8_t i; } a;
+
+  a.i = i;
+  writeOrderedBuffer(&a, 1);
+}

external/cosmotool/src/fortran.hpp

@@ -0,0 +1,113 @@
+#ifndef __COSMO_FORTRAN_HPP
+#define __COSMO_FORTRAN_HPP
+
+#include <string>
+#include <stdint.h>
+#include <inttypes.h>
+#include <iostream>
+#include "config.hpp"
+
+namespace CosmoTool
+{
+  class InvalidUnformattedAccess : public Exception
+  {
+  };
+
+  class FortranTypes
+  {
+  public:
+    enum Ordering {
+      LittleEndian, BigEndian
+    };
+
+    enum CheckpointSize {
+      Check_32bits, Check_64bits
+    };
+  };
+
+  class UnformattedRead: public FortranTypes
+  {
+  public:
+
+    UnformattedRead(const std::string& fname)
+      throw (NoSuchFileException);
+    UnformattedRead(const char *fname)
+      throw (NoSuchFileException);
+    ~UnformattedRead();
+
+    // Todo implement primitive description
+    void setOrdering(Ordering o);
+    void setCheckpointSize(CheckpointSize cs);
+    
+    void beginCheckpoint()
+      throw (InvalidUnformattedAccess,EndOfFileException);
+    void endCheckpoint(bool autodrop = false)
+      throw (InvalidUnformattedAccess);
+    
+    double readReal64()
+      throw (InvalidUnformattedAccess);
+    float readReal32()
+      throw (InvalidUnformattedAccess);
+    int32_t readInt32()
+      throw (InvalidUnformattedAccess);
+    int64_t readInt64()
+      throw (InvalidUnformattedAccess);
+
+    void skip(int64_t off)
+      throw (InvalidUnformattedAccess);
+
+  protected:
+    bool swapOrdering;
+    CheckpointSize cSize;
+    uint64_t checkPointRef;
+    uint64_t checkPointAccum;
+    std::ifstream *f;
+
+    void readOrderedBuffer(void *buffer, int size)
+      throw (InvalidUnformattedAccess);
+  };
+
+ class UnformattedWrite: public FortranTypes
+  {
+  public:
+
+    UnformattedWrite(const std::string& fname)
+      throw (NoSuchFileException);
+    UnformattedWrite(const char *fname)
+      throw (NoSuchFileException);
+    ~UnformattedWrite();
+
+    // Todo implement primitive description
+    void setOrdering(Ordering o);
+    void setCheckpointSize(CheckpointSize cs);
+    
+    void beginCheckpoint()
+      throw (FilesystemFullException,InvalidUnformattedAccess);
+    void endCheckpoint()
+      throw (FilesystemFullException,InvalidUnformattedAccess);
+    
+    void writeReal64(double d)
+      throw (FilesystemFullException);
+    void writeReal32(float f)
+      throw (FilesystemFullException);
+    void writeInt32(int32_t i)
+      throw (FilesystemFullException);
+    void writeInt64(int64_t i)
+      throw (FilesystemFullException);
+    void writeInt8(int8_t c)
+      throw (FilesystemFullException);
+
+    void writeOrderedBuffer(void *buffer, int size)
+      throw(FilesystemFullException);
+  protected:
+    bool swapOrdering;
+    CheckpointSize cSize;
+    std::streamoff checkPointRef;
+    uint64_t checkPointAccum;
+    std::ofstream *f;
+
+  };
+
+};
+
+#endif

external/cosmotool/src/fourier/base_types.hpp

@@ -0,0 +1,102 @@
+#ifndef __BASE_FOURIER_TYPES_HPP
+#define __BASE_FOURIER_TYPES_HPP
+
+#include <string>
+#include <Eigen/Dense>
+#include <complex>
+
+namespace CosmoTool
+{
+  template<typename T>
+  class FourierMap
+  {
+  protected:
+    FourierMap() {}
+    
+  public:
+    typedef Eigen::Matrix<T, 1, Eigen::Dynamic> VecType;
+    typedef Eigen::Map<VecType, Eigen::Aligned> MapType;
+    typedef Eigen::Map<VecType const, Eigen::Aligned> ConstMapType;
+
+    virtual ~FourierMap() {}
+
+    virtual const T* data() const = 0;
+    virtual T* data()  = 0;
+
+    virtual long size() const = 0;
+
+    MapType eigen()
+    {
+      return MapType(data(), size());
+    }
+
+    ConstMapType eigen() const
+    {
+      return ConstMapType(data(), size());
+    }
+ 
+    void scale(const T& factor)
+    {
+      MapType m(data(), size());
+      m *= factor;
+    }
+
+    void scale(const FourierMap<T> *map2)
+    {
+      assert(size() == map2->size());
+      MapType m(data(), size());
+      MapType m2(map2->data(), map2->size());
+      m *= m2;
+    }
+
+    void add(const T&  factor)
+    {
+      eigen() += factor;
+    }
+
+    void add(const FourierMap<T> *map2)
+    {
+      assert(size() == map2->size());
+      MapType m(data(), size());
+      MapType m2(map2->data(), map2->size());
+
+      eigen() += map2->eigen();
+    }
+
+    virtual FourierMap<T> *copy() const
+    {
+      FourierMap<T> *m = this->mimick();
+      
+      m->eigen() = this->eigen();
+      return m;
+    }
+
+    virtual FourierMap<T> *mimick() const = 0;
+  };
+
+  template<typename T>
+  class FourierTransform
+  {
+  protected:
+    FourierTransform() {}
+  public:
+    virtual ~FourierTransform() { }
+    
+    virtual const FourierMap<std::complex<T> >& fourierSpace() const = 0;
+    virtual FourierMap<std::complex<T> >& fourierSpace() = 0;
+
+    virtual const FourierMap<T>& realSpace() const = 0;
+    virtual FourierMap<T>& realSpace() = 0;
+
+    virtual FourierTransform<T> *mimick() const = 0;
+
+    virtual void analysis() = 0;
+    virtual void synthesis() = 0;
+    virtual void analysis_conjugate() = 0;
+    virtual void synthesis_conjugate() = 0;
+  };
+
+
+};
+
+#endif

external/cosmotool/src/fourier/euclidian.hpp

@@ -0,0 +1,200 @@
+#ifndef __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
+#define __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
+
+#include <vector>
+#include <boost/shared_ptr.hpp>
+#include "base_types.hpp"
+#include "fft/fftw_calls.hpp"
+
+namespace CosmoTool
+{
+  
+  template<typename T>
+  class EuclidianFourierMap: public FourierMap<T>
+  {
+  private:
+    boost::shared_ptr<T> m_data;
+    std::vector<int> m_dims;
+    long m_size;
+  public:
+    EuclidianFourierMap(boost::shared_ptr<T> indata, std::vector<int> indims)
+    {
+      m_data = indata;
+      m_dims = indims;
+      m_size = 1;
+      for (int i = 0; i < m_dims.size(); i++)
+	m_size *= m_dims[i];
+    }
+
+    virtual ~EuclidianFourierMap()
+    {
+    }
+
+    virtual const T *data() const { return m_data.get(); }
+    virtual T *data() { return m_data.get(); }
+    virtual long size() const  { return m_size; } 
+    
+    virtual FourierMap<T> *copy() const
+    {
+      FourierMap<T> *m = this->mimick();
+      m->eigen() = this->eigen();
+      return m;
+    }
+
+    virtual FourierMap<T> *mimick() const
+    {
+      return new EuclidianFourierMap<T>(
+		      boost::shared_ptr<T>((T *)fftw_malloc(sizeof(T)*m_size), 
+					   std::ptr_fun(fftw_free)),
+		       m_dims);
+    }
+  };
+
+  
+  template<typename T>
+  class EuclidianFourierTransform: public FourierTransform<T>
+  {
+  private:
+    typedef FFTW_Calls<T> calls;
+    EuclidianFourierMap<T> *realMap;
+    EuclidianFourierMap<std::complex<T> > *fourierMap;
+    typename calls::plan_type m_analysis, m_synthesis;
+    double volume;
+    long N;
+    std::vector<int> m_dims;
+    std::vector<double> m_L;
+  public:
+    EuclidianFourierTransform(const std::vector<int>& dims, const std::vector<double>& L)
+    {
+      assert(L.size() == dims.size());
+
+      m_dims = dims;
+      m_L = L;
+
+      N = 1;
+      volume = 1;
+      for (int i = 0; i < dims.size(); i++)
+	{
+	  N *= dims[i];
+	  volume *= L[i];
+	}
+
+      realMap = new EuclidianFourierMap<T>(
+		   boost::shared_ptr<T>(calls::alloc_real(N),
+					std::ptr_fun(calls::free)),
+		   dims);
+      fourierMap = new EuclidianFourierMap<std::complex<T> >(
+		   boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(N),
+						       std::ptr_fun(calls::free)), 
+		   dims);
+      m_analysis = calls::plan_dft_r2c(dims.size(), &dims[0],
+				       realMap->data(), (typename calls::complex_type *)fourierMap->data(),
+				       FFTW_MEASURE);
+      m_synthesis = calls::plan_dft_c2r(dims.size(), &dims[0],
+					(typename calls::complex_type *)fourierMap->data(), realMap->data(),
+					FFTW_MEASURE);
+    }
+    
+    virtual ~EuclidianFourierTransform()
+    {
+      delete realMap;
+      delete fourierMap;
+      calls::destroy_plan(m_synthesis);
+      calls::destroy_plan(m_analysis);
+    }
+
+    void synthesis()
+    {
+      calls::execute(m_synthesis);
+      realMap->scale(1/volume);
+    }
+    
+    void analysis()
+    {
+      calls::execute(m_analysis);
+      fourierMap->scale(volume/N);
+    }
+    
+    void synthesis_conjugate()
+    {
+      calls::execute(m_analysis);
+      fourierMap->scale(1/volume);
+    }
+    
+    void analysis_conjugate()
+    {
+      calls::execute(m_synthesis);
+      realMap->scale(volume/N);
+    }
+
+    const FourierMap<std::complex<T> >& fourierSpace() const
+    {
+      return *fourierMap;
+    }
+
+    FourierMap<std::complex<T> >& fourierSpace()
+    {
+      return *fourierMap;
+    }
+
+    const FourierMap<T>& realSpace() const
+    {
+      return *realMap;
+    }
+
+    FourierMap<T>& realSpace()
+    {
+      return *realMap;
+    }
+
+    FourierTransform<T> *mimick() const
+    {
+      return new EuclidianFourierTransform(m_dims, m_L);
+    }
+  };
+
+  template<typename T>
+  class EuclidianFourierTransform_2d: public EuclidianFourierTransform<T>
+  {
+  private:
+    template<typename T2>
+    static std::vector<T2> make_2d_vector(T2 a, T2 b)
+    {
+      T2 arr[2] = { a, b};
+      return std::vector<T2>(&arr[0], &arr[2]);
+    }
+  public:
+    EuclidianFourierTransform_2d(int Nx, int Ny, double Lx, double Ly)
+      : EuclidianFourierTransform<T>(make_2d_vector<int>(Nx, Ny), make_2d_vector<double>(Lx, Ly))
+    {
+    }
+
+    virtual ~EuclidianFourierTransform_2d() {}
+
+  };
+
+  template<typename T>
+  class EuclidianFourierTransform_3d: public EuclidianFourierTransform<T>
+  {
+  private:
+    template<typename T2>
+    static std::vector<T2> make_3d_vector(T2 a, T2 b, T2 c)
+    {
+      T2 arr[2] = { a, b, c};
+      return std::vector<T2>(&arr[0], &arr[3]);
+    }
+
+  public:
+    EuclidianFourierTransform_3d(int Nx, int Ny, int Nz, double Lx, double Ly, double Lz)
+      : EuclidianFourierTransform<T>(make_3d_vector<int>(Nx, Ny, Nz), make_3d_vector<double>(Lx, Ly, Lz))
+    {
+    }
+
+    virtual ~EuclidianFourierTransform_3d() {}
+  };
+
+
+
+};
+
+#endif

external/cosmotool/src/fourier/fft/fftw_calls.hpp

@@ -0,0 +1,75 @@
+#ifndef __FFTW_UNIFIED_CALLS_HPP
+#define __FFTW_UNIFIED_CALLS_HPP
+
+#include <fftw3.h>
+
+namespace CosmoTool
+{
+
+static inline void init_fftw_wisdom()
+{
+  fftw_import_system_wisdom();
+  fftw_import_wisdom_from_filename("fft_wisdom");
+}
+
+static inline void save_fftw_wisdom()
+{
+  fftw_export_wisdom_to_filename("fft_wisdom");
+}
+
+template<typename T> class FFTW_Calls {};
+
+
+#define FFTW_CALLS_BASE(rtype, prefix) \
+  template<>				\
+class FFTW_Calls<rtype> {		\
+public: \
+  typedef rtype real_type; \
+  typedef prefix ## _complex complex_type; \
+  typedef prefix ## _plan plan_type; \
+  \
+  static complex_type *alloc_complex(int N) { return prefix ## _alloc_complex(N); } \
+  static real_type *alloc_real(int N) { return prefix ## _alloc_real(N); } \
+  static void free(void *p) { fftw_free(p); } \
+\
+  static void execute(plan_type p) { prefix ## _execute(p); } \
+  static plan_type plan_dft_r2c_2d(int Nx, int Ny,  \
+		       real_type *in, complex_type *out, \
+		       unsigned flags) \
+  { \
+    return prefix ## _plan_dft_r2c_2d(Nx, Ny, in, out, \
+				flags); \
+  } \
+  static plan_type plan_dft_c2r_2d(int Nx, int Ny,  \
+                       complex_type *in, real_type *out, \
+                       unsigned flags) \
+  { \
+    return prefix ## _plan_dft_c2r_2d(Nx, Ny, in, out, \
+                                flags); \
+  } \
+  static plan_type plan_dft_r2c_3d(int Nx, int Ny, int Nz, \
+                                   real_type *in, complex_type *out, \
+                                   unsigned flags) \
+  { \
+    return prefix ## _plan_dft_r2c_3d(Nx, Ny, Nz, in, out, flags); \
+  } \
+  static plan_type plan_dft_r2c(int rank, const int *n, real_type *in, \
+                                complex_type *out, unsigned flags) \
+  { \
+    return prefix ## _plan_dft_r2c(rank, n, in, out, flags); \
+  } \
+  static plan_type plan_dft_c2r(int rank, const int *n, complex_type *in, \
+                                real_type *out, unsigned flags) \
+  { \
+    return prefix ## _plan_dft_c2r(rank, n, in, out, flags); \
+  } \
+  static void destroy_plan(plan_type plan) { prefix ## _destroy_plan(plan); } \
+}
+
+
+FFTW_CALLS_BASE(double, fftw);
+FFTW_CALLS_BASE(float, fftwf);
+
+};
+
+#endif

external/cosmotool/src/fourier/healpix.hpp

@@ -0,0 +1,135 @@
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_HPP
+#define __COSMOTOOL_FOURIER_HEALPIX_HPP
+
+#include <boost/bind.hpp>
+#include <boost/shared_ptr.hpp>
+#include <alm.h>
+#include <healpix_base.h>
+#include <psht_cxx.h>
+
+namespace CosmoTool
+{
+    
+  template<typename T>
+  class HealpixFourierMap: public FourierMap<T>, public Healpix_Base
+  {
+  private:
+    T *m_data;
+    Eigen::aligned_allocator<T> alloc;
+  public:
+    HealpixFourierMap(long nSide)
+      : Healpix_Base(RING, nSide, SET_NSIDE)
+    {
+      m_data = alloc.allocate(Npix);
+    }
+
+    virtual ~HealpixFourierMap()
+    {
+      alloc.deallocate(m_data);
+    }
+
+    virtual const T* data() const { return m_data; }
+    virtual T *data() { return m_data; }
+    virtual long size() const { return Npix(); }
+    
+    virtual FourierMap<T> *mimick() const
+    {
+      return new HealpixFourierMap<T>(Nside());
+    }
+  };
+
+
+  template<typename T>
+  class HealpixFourierALM: public FourierMap<std::complex<T> >, public Alm_Base
+  {
+  private:
+    std::complex<T> *alms;
+    long size;
+    Eigen::aligned_allocator<std::complex<T> > alloc;
+  public:
+    HealpixFourierALM(long Lmax, long Mmax)
+      : Alm_Base(Lmax, Mmax)
+    {
+      size = Num_Alms(Lmax, Mmax);
+      alms = alloc.allocate(size);
+    }
+
+    virtual ~HealpixFourierALM()
+    {
+      alloc.deallocate(alms);
+    }
+
+    virtual const T* data() const { return alms; }
+    virtual T * data() { return alms;} 
+    virtual long size() const { return size; }
+
+    virtual FourierMap<T> *mimick() const
+    {
+      return new HealpixFourierALM<T>(Lmax(), Mmax());
+    }
+  };
+
+  template<typename T>
+  class HealpixFourierTransform: public FourierTransform<T>
+  {
+  private:
+    HealpixFourierMap<T> realMap;
+    HealpixFourierALM<T> fourierMap;
+    psht_joblist<T> jobs;
+  public:
+    HealpixFourierTransform(long nSide, long Lmax, long Mmax)
+      : realMap(nSide), fourierMap(Lmax, Mmax)
+    {
+      jobs.set_Healpix_geometry(nSide);
+      jobs.set_triangular_alm_info(Lmax, Mmax);
+    }
+
+    virtual ~HealpixFourierTransform() {}
+
+    virtual const FourierMap<std::complex<T> >& fourierSpace() const { return fourierMap; }
+
+    virtual FourierMap<std::complex<T> >& fourierSpace() { return fourierMap; }
+
+    virtual const FourierMap<T>& realSpace() const { return realMap; }
+
+    virtual FourierMap<T>& realSpace() { return realMap; }
+    
+    virtual FourierTransform<T> *mimick() const 
+    {
+      return new HealpixFourierTransform<T>(realMap.Nside(), fourierMap.Lmax(), fourierMap.Mmax());
+    }
+
+    virtual void analysis()
+    {
+      jobs.add_map2alm(realMap.data(), 
+		       reinterpret_cast<xcomplex<T> *>(fourierMap.data()), 
+		       false);
+      jobs.execute();
+      jobs.clear_jobs();
+    }
+
+    virtual void synthesis()
+    {
+      jobs.add_alm2map(reinterpret_cast<xcomplex<T> *>(fourierMap.data()), 
+		       realMap.data(),
+		       false);
+      jobs.execute();
+      jobs.clear_jobs();
+    }
+
+    virtual void analysis_conjugate()
+    {
+      synthesis();
+      realMap.scale(4*M_PI/realMap.Npix());
+    }
+
+    virtual void synthesis_conjugate()
+    {
+      analysis();
+      fourierMap.scale(realMap.Npix()/(4*M_PI));
+    }
+    
+  };
+};
+
+#endif

external/cosmotool/src/growthFactor.cpp

@@ -0,0 +1,111 @@
+#include <cmath>
+#include <gsl/gsl_integration.h>
+#include "interpolate.hpp"
+#include "growthFactor.hpp"
+
+using namespace CosmoTool;
+
+#define AMIN 1e-5
+#define AMAX 1.0
+#define NUM_WORK 5000
+#define TOLERANCE 1e-6
+
+typedef struct {
+  double OmegaLambda;
+  double OmegaMatter;
+  double Hubble;
+} Cosmology;
+
+static double computeOmegaMatter(Cosmology *cosmo, double a)
+{
+  return cosmo->OmegaMatter / (cosmo->OmegaMatter + a*a*a * cosmo->OmegaLambda);
+}
+
+static double computeHdotH(Cosmology *cosmo, double a)
+{
+  return -1.5 * cosmo->OmegaMatter / (a * (cosmo->OmegaMatter + a*a*a*cosmo->OmegaLambda));
+}
+
+static double computeE(double OmegaMatter, double OmegaLambda, double a)
+{
+  double H2;
+  double OmegaK = (1 - OmegaMatter - OmegaLambda);
+  
+  H2 = OmegaMatter/(a*a*a) + OmegaLambda + OmegaK/(a*a);
+
+  return sqrt(H2);
+}
+
+static double computeEprime(Cosmology *cosmo, double a)
+{
+  double H2;
+  double OmegaK = (1 - cosmo->OmegaMatter - cosmo->OmegaLambda);
+  
+  H2 = -3*cosmo->OmegaMatter/(a*a*a*a) - 2*OmegaK/(a*a*a);
+
+  return 0.5*H2/computeE(cosmo->OmegaMatter, cosmo->OmegaLambda, a);  
+}
+
+static inline double cube(double x)
+{
+  return x*x*x;
+}
+
+static double integrandGrowthFactor(double a, void *params)
+{
+  Cosmology *cosmo = (Cosmology *)params;
+
+  return 1/cube(computeE(cosmo->OmegaMatter, cosmo->OmegaLambda, a)*a);
+}
+
+Interpolate CosmoTool::buildLinearGrowth(double OmegaLambda, double OmegaMatter, double Hubble, int numPoints)
+{
+  Cosmology cosmology;
+  gsl_integration_workspace *work = gsl_integration_workspace_alloc(NUM_WORK);
+  gsl_function f;
+  double *a_input, *D_output;
+
+  cosmology.OmegaLambda = OmegaLambda;
+  cosmology.OmegaMatter = OmegaMatter;
+  cosmology.Hubble = Hubble;
+
+  a_input = new double[numPoints];
+  D_output = new double[numPoints];
+
+  f.params = &cosmology;
+  f.function = integrandGrowthFactor;
+  
+  a_input[0] = 0;
+  D_output[0] = 0;
+
+  for (int i = 1; i < numPoints; i++)
+    {
+      double a_dest = 0 + 1.0*i/(numPoints-1);
+      double result, abserr;
+      double E = computeE(cosmology.OmegaMatter, cosmology.OmegaLambda, a_dest);
+      double Eprime = computeEprime(&cosmology, a_dest);
+
+      gsl_integration_qag(&f, 0, a_dest, 0, TOLERANCE, NUM_WORK, 
+			  GSL_INTEG_GAUSS61, work, &result, &abserr);
+
+
+      result *= 2.5 * computeE(cosmology.OmegaMatter, cosmology.OmegaLambda, a_dest) * OmegaMatter;
+
+      D_output[i] = result;
+      a_input[i] = a_dest;
+
+    }
+  gsl_integration_workspace_free(work);
+
+  for (int i = 0; i < numPoints; i++)
+    {
+      D_output[i] /= D_output[numPoints-1];
+    }
+
+  Interpolate p(a_input, D_output, numPoints, true, false, true);
+
+  delete[] a_input;
+  delete[] D_output;
+
+  return p;
+}

external/cosmotool/src/growthFactor.hpp

@@ -0,0 +1,12 @@
+#ifndef COSMO_GROWTH_FACTOR_HPP
+#define COSMO_GROWTH_FACTOR_HPP
+
+#include "interpolate.hpp"
+
+namespace CosmoTool
+{
+  Interpolate buildLinearGrowth(double OmegaLambda, double OmegaMatter, double Hubble, int numPoints = 10000);
+
+};
+
+#endif

external/cosmotool/src/h5_readFlash.cpp

@@ -0,0 +1,474 @@
+/* This file contains the functions that read the data from the HDF5 file
+ * The functions accept the PARAMESH data through arguments.
+ */
+
+#include "h5_readFlash.hpp"
+
+using namespace H5;
+
+/* indices of attributes in fof_particle_type */
+int iptag_out = 0;
+int ipx_out   = 0;
+int ipy_out   = 1;
+int ipz_out   = 2;
+int ipvx_out   = 3;
+int ipvy_out   = 4;
+int ipvz_out   = 5;
+
+/* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx */
+
+/*  n*_runtime_parameters should be set by the caller to 
+    the maximum number of runtime parameters to read. 
+*/
+
+void h5_read_runtime_parameters
+     (H5File* file,    /* file handle */
+      double* LBox,    // box size
+      int* numPart,
+      double *hubble,
+      double *omegam,
+      double *omegalambda)       // number of particles
+{
+
+ int MAX_PARM = 200;
+
+ int nreal_runtime_parameters, nint_runtime_parameters, nstr_runtime_parameters;
+
+ char real_runtime_parameter_names[MAX_PARM][RUNTIME_PARAMETER_STRING_SIZE];
+ char int_runtime_parameter_names[MAX_PARM][RUNTIME_PARAMETER_STRING_SIZE];
+ char str_runtime_parameter_names[MAX_PARM][RUNTIME_PARAMETER_STRING_SIZE];
+
+ double real_runtime_parameter_values[MAX_PARM];
+ int    int_runtime_parameter_values[MAX_PARM];
+ char   str_runtime_parameter_values[MAX_PARM][RUNTIME_PARAMETER_STRING_SIZE];
+
+  int rank;
+  hsize_t dimens_1d, maxdimens_1d;
+
+  real_runtime_params_t *real_rt_parms;
+  int_runtime_params_t *int_rt_parms;
+  str_runtime_params_t *str_rt_parms;
+  log_runtime_params_t *log_rt_parms;
+
+  double omegarad;
+
+  int i;
+
+  nint_runtime_parameters = MAX_PARM;
+  nreal_runtime_parameters = MAX_PARM;
+
+  /* integer runtime parameters */
+  int_rt_parms = (int_runtime_params_t *) malloc(nint_runtime_parameters * sizeof(int_runtime_params_t)); 
+
+  rank = 1;
+  DataSet dataset = file->openDataSet("integer runtime parameters");
+  
+  IntType int_rt_type = dataset.getIntType();  
+  //int_rt_type = H5Dget_type(dataset);
+
+  DataSpace dataspace = dataset.getSpace();
+  //dataspace = H5Dget_space(dataset);
+
+  int ndims = dataspace.getSimpleExtentDims(&dimens_1d, NULL);
+  //H5Sget_simple_extent_dims(dataspace, &dimens_1d, &maxdimens_1d);
+
+  /* don't read in more than we can handle */
+  if (nint_runtime_parameters < dimens_1d) {
+    dimens_1d = nint_runtime_parameters;
+  } else {
+    nint_runtime_parameters = dimens_1d;
+  }
+  DataSpace memspace(rank, &dimens_1d);
+  //memspace = H5Screate_simple(rank, &dimens_1d, NULL);
+
+  dataset.read(int_rt_parms, int_rt_type, memspace, dataspace,
+		    H5P_DEFAULT);
+  //status = H5Dread(dataset, int_rt_type, memspace, dataspace,
+	//	    H5P_DEFAULT, int_rt_parms);
+
+
+  for (i = 0; i < nint_runtime_parameters; i++) {
+    strncpy(int_runtime_parameter_names[i], int_rt_parms[i].name, RUNTIME_PARAMETER_STRING_SIZE);
+    int_runtime_parameter_values[i] = int_rt_parms[i].value;
+  }
+
+  free(int_rt_parms);
+  memspace.close();
+  dataspace.close();
+  dataset.close();
+  //H5Sclose(dataspace);
+  //H5Dclose(dataset);
+
+  /* done with int runtime parameters */
+
+  /* reals */
+
+  real_rt_parms = (real_runtime_params_t *) malloc(nreal_runtime_parameters * sizeof(real_runtime_params_t)); 
+  
+  rank = 1;
+  dataset = file->openDataSet("real runtime parameters");
+  //dataset = H5Dopen(*file_identifier, "real runtime parameters"); 
+  
+  dataspace = dataset.getSpace();
+  FloatType real_rt_type = dataset.getFloatType();  
+  ndims = dataspace.getSimpleExtentDims(&dimens_1d, NULL);
+  //dataspace = H5Dget_space(dataset);
+  //real_rt_type = H5Dget_type(dataset);
+  //H5Sget_simple_extent_dims(dataspace, &dimens_1d, &maxdimens_1d);
+
+  /* don't read in more than we can handle */
+  if (nreal_runtime_parameters < dimens_1d) {
+    dimens_1d = nreal_runtime_parameters;
+  } else {
+    nreal_runtime_parameters = dimens_1d;
+  }
+   memspace = DataSpace(rank, &dimens_1d);
+  //memspace = H5Screate_simple(rank, &dimens_1d, NULL);
+
+  dataset.read(real_rt_parms, real_rt_type, memspace, dataspace,
+		    H5P_DEFAULT);
+  //status = H5Dread(dataset, real_rt_type, memspace, dataspace,
+//		    H5P_DEFAULT, real_rt_parms);
+
+
+  for (i = 0; i < nreal_runtime_parameters; i++) {
+    strncpy(real_runtime_parameter_names[i], real_rt_parms[i].name, RUNTIME_PARAMETER_STRING_SIZE);
+    real_runtime_parameter_values[i] = real_rt_parms[i].value;
+  }
+
+  free(real_rt_parms);
+  memspace.close();
+  dataspace.close();
+  dataset.close();
+  //H5Sclose(dataspace);
+  //H5Dclose(dataset);
+
+  /* done with reals */
+
+  // grab the data we want
+  for (i = 0; i < nreal_runtime_parameters; i++) {
+    if (strncmp(real_runtime_parameter_names[i],"xmax",4) == 0 ) {
+      *LBox = real_runtime_parameter_values[i];
+    }
+    if (strncmp(real_runtime_parameter_names[i],"hubbleconstant", 14) == 0 ) {
+      *hubble = real_runtime_parameter_values[i];      
+    }
+    if (strncmp(real_runtime_parameter_names[i],"omegamatter", 11) == 0 ) {
+      *omegam = real_runtime_parameter_values[i];
+    }
+    if (strncmp(real_runtime_parameter_names[i],"omegaradiation", 11) == 0 ) {
+      omegarad = real_runtime_parameter_values[i];
+    }
+    if (strncmp(real_runtime_parameter_names[i],"cosmologicalconstant", 20) == 0 ) {
+      *omegalambda = real_runtime_parameter_values[i];
+    }
+  }
+  
+  for (i = 0; i < nint_runtime_parameters; i++) {
+    if (strncmp(int_runtime_parameter_names[i],"pt_numx",7) == 0 ) {
+      *numPart = int_runtime_parameter_values[i];
+      *numPart = *numPart * *numPart * *numPart;
+    }
+  }
+
+}
+
+     
+/* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
+void h5_read_flash3_particles (H5File* file,
+                        int* totalparticles,
+                        int* localnp,
+                        int* particle_offset,
+                        float *pos1,
+                        float *pos2,
+                        float *pos3,
+                        float *vel1,
+                        float *vel2,
+                        float *vel3,
+                        int    *id)
+{
+
+  herr_t   status;
+  hsize_t  dimens_1d, maxdims_1d;
+  hsize_t  start_1d, stride_1d, count_1d;
+  int rank;
+
+  int      numProps, i, p;
+  int      numPartBuffer = 5000, sizePartBuffer, sizePart;
+  int      pstack, poffset, pcount;
+  int      iptag, ipx, ipy, ipz, ipvx, ipvy, ipvz;
+  char     *propName;
+  double   *partBuffer;
+  char    part_names[50][OUTPUT_PROP_LENGTH];
+  int     string_size;  
+         
+//  char part_names[NPART_PROPS][OUTPUT_PROP_LENGTH];
+  hsize_t dimens_2d[2], maxdimens_2d[2];
+  hsize_t start_2d[2], count_2d[2], stride_2d[2];
+
+  /* skip this routine if no particles to read */
+  if ((*localnp) == 0) {
+    return;
+  }
+
+
+ /* first determine how many particle properties are
+     present in the input data file, and determine which of these
+     are the properties we are interested in */
+  DataSet dataset = file->openDataSet("particle names");
+  DataSpace dataspace = dataset.getSpace();
+  //dataset = H5Dopen(*file_identifier, "particle names");
+  //dataspace = H5Dget_space(dataset);
+
+  int ndims = dataspace.getSimpleExtentDims(dimens_2d, NULL);
+  //H5Sget_simple_extent_dims(dataspace, dimens_2d, maxdimens_2d);
+
+  //total number of particle properties
+  numProps = dimens_2d[0];
+  
+  string_size = OUTPUT_PROP_LENGTH;
+  StrType string_type = H5Tcopy(H5T_C_S1);
+  string_type.setSize(string_size);
+  //status = H5Tset_size(string_type, string_size);
+  
+  rank = 2;
+
+  start_2d[0] = 0;
+  start_2d[1] = 0;
+
+  stride_2d[0] = 1;
+  stride_2d[1] = 1;
+
+  count_2d[0] = dimens_2d[0];
+  count_2d[1] = dimens_2d[1];
+
+  dataspace.selectHyperslab(H5S_SELECT_SET, count_2d, start_2d); 
+  //status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, start_2d,
+  //                      stride_2d, count_2d, NULL);
+
+  DataSpace memspace(rank, dimens_2d);
+  //memspace = H5Screate_simple(rank, dimens_2d, NULL);
+
+  dataset.read(part_names, string_type, H5S_ALL, H5S_ALL, H5P_DEFAULT);
+  //status = H5Dread(dataset, string_type, H5S_ALL, H5S_ALL,
+  //              H5P_DEFAULT, part_names);
+
+ 
+  string_type.close();
+  memspace.close();
+  dataspace.close();
+  dataset.close(); 
+  //H5Tclose(string_type);
+  //H5Sclose(memspace);
+  //H5Sclose(dataspace);
+  //H5Dclose(dataset);
+
+  for (i=0;i<numProps;i++) {
+    if (strncmp(part_names[i], "tag", 3) == 0)  { iptag = i+1; }
+    if (strncmp(part_names[i], "posx", 4) == 0) { ipx = i+1; }
+    if (strncmp(part_names[i], "posy", 4) == 0) { ipy = i+1; }
+    if (strncmp(part_names[i], "posz", 4) == 0) { ipz = i+1; }
+    if (strncmp(part_names[i], "velx", 4) == 0) { ipvx = i+1; }
+    if (strncmp(part_names[i], "vely", 4) == 0) { ipvy = i+1; }
+    if (strncmp(part_names[i], "velz", 4) == 0) { ipvz = i+1; }
+  }
+
+  if ((iptag < 0) || (ipx < 0) || (ipy < 0) || (ipz < 0) || (ipvx < 0) || 
+       (ipvy < 0) || (ipvz < 0) ) {
+    printf("One or more required particle attributes not found in file!\n");
+    return;
+  }
+
+  //printf("iptag = %d, ipx = %d, ipy = %d, ipz = %d\n", iptag, ipx, ipy, ipz);
+  //printf("ipvx = %d, ipvy = %d, ipvz = %d\n", ipvx, ipvy, ipvz);
+
+  //read particles
+  dataset = file->openDataSet("tracer particles");
+  //dataset = H5Dopen(*file_identifier, "tracer particles");
+  
+  FloatType datatype = dataset.getFloatType();
+  //datatype = H5Dget_type(dataset);
+
+  /* establish read-in particle buffer */
+  sizePart = numProps*(sizeof(double));
+  sizePartBuffer = numPartBuffer * sizePart;
+  partBuffer = (double *)malloc(sizePartBuffer);
+
+  dataspace = dataset.getSpace();
+  ndims = dataspace.getSimpleExtentDims(dimens_2d, NULL);
+  //dataspace = H5Dget_space(dataset);
+  //H5Sget_simple_extent_dims(dataspace, dimens_2d, maxdimens_2d);
+ 
+  /*insert particle properties (numPartBuffer) particles at a time*/
+  pstack = (*localnp);
+  poffset = 0;
+  if (pstack > numPartBuffer) {
+   pcount = numPartBuffer;
+  }
+  else {
+   pcount = pstack;
+  }
+
+  while ( pstack > 0) {
+    rank       = 2;
+    maxdimens_2d[0] = (hsize_t) (*totalparticles);
+    maxdimens_2d[1] = (hsize_t) (numProps);
+
+    start_2d[0]  = (hsize_t) (*particle_offset + poffset);
+    start_2d[1]  = (hsize_t) 0;
+
+    stride_2d[0] = 1;
+    stride_2d[1] = 1;
+
+    count_2d[0]  = (hsize_t) (pcount);
+    count_2d[1]  = (hsize_t) (numProps);
+
+    dimens_2d[0] = (pcount);
+    dimens_2d[1] = (numProps);
+
+    dataspace.selectHyperslab(H5S_SELECT_SET, count_2d, start_2d); 
+    //status     = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, start_2d,
+    //                                 stride_2d, count_2d, NULL);
+
+    memspace = DataSpace(rank, dimens_2d);
+    //memspace   = H5Screate_simple(rank, dimens_2d, maxdimens_2d);
+
+    /* read data from the dataset */
+   dataset.read(partBuffer, datatype, memspace, dataspace, H5P_DEFAULT);
+   //status = H5Dread(dataset, datatype, memspace, dataspace, H5P_DEFAULT, partBuffer);
+
+    /* convert buffer into particle struct */
+
+    if (id) {
+    for(p=0; p < (pcount); p++) {
+      id[p+poffset] = (int)  *(partBuffer+iptag-1+p*numProps);
+    } }
+
+    if (pos1 && pos2 && pos3) {
+    for(p=0; p < (pcount); p++) {
+      pos1[p+poffset]  = (float) *(partBuffer+ipx-1+p*numProps);
+      pos2[p+poffset]  = (float) *(partBuffer+ipy-1+p*numProps);
+      pos3[p+poffset]  = (float) *(partBuffer+ipz-1+p*numProps);
+    }
+    }
+   
+    
+    if (vel1 && vel2 && vel3) {
+    for(p=0; p < (pcount); p++) {
+      vel1[p+poffset]  = (float) *(partBuffer+ipvx-1+p*numProps);
+      vel2[p+poffset]  = (float) *(partBuffer+ipvy-1+p*numProps);
+      vel3[p+poffset]  = (float) *(partBuffer+ipvz-1+p*numProps);
+    }
+    }
+
+    memspace.close();
+    //status = H5Sclose(memspace);
+   /* advance buffer */
+    pstack = pstack - pcount;
+    poffset = poffset + pcount;
+    if (pstack > numPartBuffer) {
+      pcount = numPartBuffer;
+    }
+    else {
+      pcount = pstack;
+    }
+
+  } /* end while */
+
+
+  datatype.close();
+  dataspace.close();
+  dataset.close();
+  //status = H5Tclose(datatype);
+  //status = H5Sclose(dataspace);
+  //status = H5Dclose(dataset);
+
+}
+ 
+
+/*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
+
+void h5_read_flash3_header_info(H5File* file,
+				double* time,                  /* simulation time */
+				double* redshift)              /* redshift of checkpoint */
+{
+
+  herr_t status;
+
+  int file_version;
+
+  hid_t sp_type, si_type;
+
+  hsize_t dimens_1d, maxdimens_1d;
+  hid_t string_type;
+  real_list_t *real_list;
+  int* num_real, num_int;
+  int MAX_SCALARS = 100;
+  char real_names[MAX_SCALARS][MAX_STRING_LENGTH];
+  double real_values[MAX_SCALARS];
+  char int_names[MAX_SCALARS][MAX_STRING_LENGTH];
+  int int_values[MAX_SCALARS];
+  int i;
+
+  H5std_string DATASET_NAME;
+
+  string_type = H5Tcopy(H5T_C_S1);  
+  H5Tset_size(string_type, MAX_STRING_LENGTH);
+
+  DataSet dataset = file->openDataSet("real scalars");
+  DataSpace dataspace = dataset.getSpace();
+
+  /* read extent of 'dataspace' (i.e. # of name/value pairs) into 'dimens_1d' */
+  int ndims = dataspace.getSimpleExtentDims(&dimens_1d, NULL);
+
+  if (dimens_1d > MAX_SCALARS) {
+    printf("Error: reading more than MAX_SCALARS runtime parameters in checkpoint file!\n");
+  }
+
+  /* malloc a pointer to a list of real_list_t's */
+  real_list = (real_list_t *) malloc(dimens_1d * sizeof(real_list_t));
+
+  // create a new simple dataspace of 1 dimension and size of 'dimens_1d' 
+  DataSpace memspace(1, &dimens_1d);
+
+  // create an empty vessel sized to hold one real_list_t's worth of data 
+  CompType real_list_type( sizeof(real_list_t) );
+
+  // subdivide the empty vessel into its component sections (name and value) 
+  real_list_type.insertMember(
+          "name",
+          HOFFSET(real_list_t, name),
+          string_type);
+
+  real_list_type.insertMember(
+          "value",
+          HOFFSET(real_list_t, value),
+          PredType::NATIVE_DOUBLE);
+
+  // read the data into 'real_list' 
+  dataset.read( real_list, real_list_type, memspace, dataspace,
+                H5P_DEFAULT);
+
+
+  if (status < 0) {
+    printf("Error readingruntime parameterss from data file\n");
+  }
+
+  for (i = 0; i < dimens_1d; i++) {
+    strncpy(real_names[i], real_list[i].name, MAX_STRING_LENGTH);
+    real_values[i] = real_list[i].value;
+
+    if (strncmp(real_names[i],"time",4) == 0 ) {
+      *time = real_values[i];
+    }
+    if (strncmp(real_names[i],"redshift",8) == 0 ) {
+      *redshift = real_values[i];
+    }
+  }
+
+  free(real_list);
+  real_list_type.close();
+  memspace.close();
+  dataspace.close();
+  dataset.close();
+}

external/cosmotool/src/h5_readFlash.hpp

@@ -0,0 +1,36 @@
+/* This file contains the functions that read the data from the HDF5 file
+ * The functions accept the PARAMESH data through arguments.
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stddef.h>
+#include <string.h>
+#include "hdf5_flash.h"
+#include "H5Cpp.h"
+
+using namespace H5;
+
+void h5_read_runtime_parameters
+     (H5File* file,    /* file handle */
+      double* LBox,
+      int* numPart,
+      double* hubble,
+      double* omegam,
+      double* omegalambda);
+
+void h5_read_flash3_particles (H5File* file,
+                        int* totalparticles,
+                        int* localnp,
+                        int* particle_offset,
+                        float *pos1,
+                        float *pos2,
+                        float *pos3,
+                        float *vel1,
+                        float *vel2,
+                        float *vel3,
+                        int    *id);
+
+void h5_read_flash3_header_info(H5File* file,
+				double* time,                   /* simulation time */
+				double *redshift);              /* simulation redshift */

external/cosmotool/src/hdf5_flash.h

@@ -0,0 +1,192 @@
+/* general header file for the HDF 5 IO in FLASH */
+
+
+#ifndef _HDF5_FLASH_H
+#define _HDF5_FLASH_H
+
+/* pull in some basic FLASH information */
+
+//#include "flash_defines.fh"
+
+/* define an integer file format version number, that is stored
+   in the output files.  This way, people can check this number
+   before reading, and compare against a published format to know
+   what is stored in the file.  In theory, this number should be
+   incremented anytime a change is made to the file format */
+
+/* File format history:
+
+   1 -- original version
+
+   2 -- added build records:
+         "FLASH build date"
+         "FLASH build directory"
+         "FLASH build machine"
+         "FLASH setup call"
+
+        added the "file format version" record
+
+   3 -- added the "run comment record" (why was this
+        not done long ago?)
+
+   4 -- added extrema attributes to the variable records
+
+   5 -- redshift included
+
+   6 -- added the Module data to the attributes of "/"
+
+   7 -- make build info attributes on "/"
+*/
+
+#define FILE_FORMAT_VERSION 7
+
+#define RUNTIME_PARAMETER_STRING_SIZE 80
+
+#define TIMER_NAME_STRING_SIZE 30
+
+#define MAX_STRING_LENGTH 80
+
+#define LIST_STRING_SIZE 80
+
+#define OUTPUT_PROP_LENGTH 24
+
+typedef struct real_list_t {
+  char name[LIST_STRING_SIZE];
+  double value;
+} real_list_t;
+
+typedef struct int_runtime_params_t {
+  char name[RUNTIME_PARAMETER_STRING_SIZE];
+  int value;
+} int_runtime_params_t;
+
+typedef struct real_runtime_params_t {
+  char name[RUNTIME_PARAMETER_STRING_SIZE];
+  double value;
+} real_runtime_params_t;
+
+typedef struct str_runtime_params_t {
+  char value[RUNTIME_PARAMETER_STRING_SIZE];
+  char name[RUNTIME_PARAMETER_STRING_SIZE];
+} str_runtime_params_t;
+
+typedef struct log_runtime_params_t {
+  int value;
+  char name[RUNTIME_PARAMETER_STRING_SIZE];
+} log_runtime_params_t;
+
+
+#define MAX_TIMER_PARENTS 20
+#define MAX_TIMER_CALL_STACK_DEPTH 20
+typedef struct timer_data_t {
+  char name[TIMER_NAME_STRING_SIZE];
+  double t_value[MAX_TIMER_PARENTS];
+  int    t_counts[MAX_TIMER_PARENTS];
+  int    t_on[MAX_TIMER_PARENTS];
+  int    t_stacks[MAX_TIMER_PARENTS][MAX_TIMER_CALL_STACK_DEPTH];
+  int    t_num_parents;
+  int    t_stack_sizes[MAX_TIMER_PARENTS];
+} full_timer_data_t;
+
+
+typedef struct sim_params_t {
+  int total_blocks;
+  int nsteps;
+  int nxb;
+  int nyb;
+  int nzb;
+  double time; 
+  double timestep;
+  double redshift;
+    
+} sim_params_t;
+
+typedef struct sim_params_sp_t {
+  int total_blocks;
+  int nsteps;
+  int nxb;
+  int nyb;
+  int nzb;
+  float time; 
+  float timestep;
+  float redshift;
+    
+} sim_params_sp_t;
+
+typedef struct sim_info_t {
+  int file_format_version;
+  char setup_call[400];
+  char file_creation_time[MAX_STRING_LENGTH];
+  char flash_version[MAX_STRING_LENGTH];
+  char build_date[MAX_STRING_LENGTH];
+  char build_dir[MAX_STRING_LENGTH];
+  char build_machine[MAX_STRING_LENGTH];
+  char cflags[400];
+  char fflags[400];
+  char setup_time_stamp[MAX_STRING_LENGTH];
+  char build_time_stamp[MAX_STRING_LENGTH];
+} sim_info_t;
+
+/* define some particle property constants */
+
+#if FLASH_NUMBER_OF_INT_PARTICLE_PROPS > 0
+#define NUMINTPROPS  2*((FLASH_NUMBER_OF_INT_PARTICLE_PROPS+1)/2)
+#else
+#define NUMINTPROPS  2
+#endif
+
+#if FLASH_NUMBER_OF_REAL_PARTICLE_PROPS > 0
+#define NUMREALPROPS FLASH_NUMBER_OF_REAL_PARTICLE_PROPS
+#else
+#define NUMREALPROPS 1
+#endif
+
+#define NUMACTUALINTPROPS FLASH_NUMBER_OF_INT_PARTICLE_PROPS
+#define NUMACTUALREALPROPS FLASH_NUMBER_OF_REAL_PARTICLE_PROPS
+
+
+/* set the dimension and grid variables -- the variable N_DIM is set 
+   in the compile line */
+
+
+/* mdim is the maximum dimension -- this is set in tree.fh */
+#define MDIM 3
+#define MGID 15
+
+/* 3-d problem */
+#if N_DIM == 3 
+
+#define NDIM  3
+
+#define NGID 15
+
+#define k2d 1
+#define k3d 1
+
+
+/* 2-d problem */
+#elif N_DIM == 2
+
+#define NDIM  2
+
+#define NGID 9
+
+#define k2d 1
+#define k3d 0
+
+
+/* 1-d problem */
+#else
+
+#define NDIM 1
+
+#define NGID 5
+
+#define k2d 0
+#define k3d 0
+
+#endif
+
+
+#endif
+

external/cosmotool/src/interpolate.cpp

@@ -0,0 +1,278 @@
+#include <cassert>
+#include <cmath>
+#include <algorithm>
+#include "interpolate.hpp"
+#include <vector>
+#include <fstream>
+#include <inttypes.h>
+#include <sstream>
+
+using namespace std;
+using namespace CosmoTool;
+  
+Interpolate::Interpolate(double *xs, double *values, uint32_t N, bool autofree,
+			 bool logx, bool logy)
+{
+  spline = gsl_spline_alloc (gsl_interp_linear, N);
+  gsl_spline_init (spline, xs, values, N);
+  accel_interp = gsl_interp_accel_alloc();
+
+  this->logx = logx;
+  this->logy = logy;
+  this->autoFree = autofree;
+}
+
+Interpolate::~Interpolate()
+{
+  if (spline != 0)
+    gsl_spline_free (spline);
+  if (accel_interp != 0)
+    gsl_interp_accel_free (accel_interp);
+}
+
+double Interpolate::compute(double x)
+  throw (InvalidRangeException)
+{
+  double y;
+  
+  if (logx)
+    x = log(x);
+
+  int err = gsl_spline_eval_e(spline, x, accel_interp, &y);
+
+  if (err)
+    throw InvalidRangeException("Interpolate argument outside range");
+
+  if (logy)
+    return exp(y);
+  else
+    return y;
+}
+
+double Interpolate::compute(double x) const
+  throw (InvalidRangeException)
+{
+  double y;
+
+  if (logx)
+    x = log(x);
+
+  int err = gsl_spline_eval_e(spline, x, 0, &y);
+
+  if (err)
+    throw InvalidRangeException("Interpolate argument outside range");
+
+  if (logy)
+    return exp(y);
+  else
+    return y;
+}
+
+
+
+
+double Interpolate::derivative(double x)
+  throw (InvalidRangeException)
+{
+  double y, dy, x0 = x;
+  
+  if (logx)
+    x0 = log(x0);
+
+  int err = gsl_spline_eval_deriv_e(spline, x0, accel_interp, &dy);
+
+  if (err)
+    throw InvalidRangeException("Interpolate argument outside range");
+
+  if (logy)
+    {
+      int err = gsl_spline_eval_e(spline, x0, accel_interp, &y);
+      
+      assert(err == 0);
+      
+      return dy*exp(y)/x0;
+    }
+  else
+    return dy;
+}
+
+uint32_t Interpolate::getNumPoints() const
+{
+  return spline->size;
+}
+
+void Interpolate::fillWithXY(double *x, double *y) const
+{
+  if (x != 0)
+    memcpy(x, spline->x, sizeof(double)*spline->size);
+  if (y != 0)
+    memcpy(y, spline->y, sizeof(double)*spline->size);
+}
+
+const Interpolate& Interpolate::operator=(const Interpolate& a)
+{
+  double *x, *y;
+
+  if (spline != NULL)
+    {
+      gsl_spline_free (spline);
+      gsl_interp_accel_free (accel_interp);
+    }
+
+  autoFree = true;
+  spline = gsl_spline_alloc (gsl_interp_linear, a.spline->size);
+  accel_interp = gsl_interp_accel_alloc ();
+  gsl_spline_init(spline, a.spline->x, a.spline->y, a.spline->size);
+  logx = a.logx;
+  logy = a.logy;
+}
+
+double Interpolate::getMaxX() const
+{
+  if (logx)
+    return exp(spline->x[spline->size-1]);
+  else
+    return spline->x[spline->size-1];
+}
+
+typedef struct {
+  double a, b;
+} MyPair;
+
+bool operator<(const MyPair& a, const MyPair& b)
+{
+  return a.a < b.a;
+}
+
+Interpolate CosmoTool::buildFromVector(const InterpolatePairs& v)
+{
+  double *x = new double[v.size()];
+  double *y = new double[v.size()];
+  
+  for (uint32_t i = 0; i < v.size(); i++)
+    {
+      x[i] = v[i].first;
+      y[i] = v[i].second;
+    }
+
+  Interpolate inter = Interpolate(x, y, v.size(), true);
+
+  delete[] x;
+  delete[] y;
+
+  return inter;  
+}
+
+Interpolate CosmoTool::buildInterpolateFromFile(const char *fname)
+  throw (NoSuchFileException)
+{
+  vector<MyPair> allData;
+  ifstream f(fname);
+
+  if (!f)
+    throw NoSuchFileException(fname);
+
+  do 
+    {
+      MyPair m;
+
+      if (!(f >> m.a >> m.b))
+	break;
+
+      allData.push_back(m);
+    }
+  while (1);
+
+  sort(allData.begin(), allData.end());
+
+  double *x = new double[allData.size()];
+  double *y = new double[allData.size()];
+  
+  for (uint32_t i = 0; i < allData.size(); i++)
+    {
+      x[i] = allData[i].a;
+      y[i] = allData[i].b;
+    }
+
+  Interpolate inter = Interpolate(x, y, allData.size(), true);
+
+  delete[] x;
+  delete[] y;
+
+  return inter;
+}
+
+Interpolate CosmoTool::buildInterpolateFromColumns(const char *fname, uint32_t col1, uint32_t col2, bool logx,
+						   bool logy)
+  throw (NoSuchFileException,InvalidRangeException)
+{
+  vector<MyPair> allData;
+  ifstream f(fname);
+
+  if (!f)
+    throw NoSuchFileException(fname);
+  
+  bool swapped = (col1 > col2);
+  uint32_t colMin = min(col1, col2);
+  uint32_t colMax = max(col1, col2);
+
+  do 
+    {
+      MyPair m;
+      string line;
+
+      if (getline(f, line).eof())
+	break;
+
+      istringstream iss(line);
+      double dummy;
+      double val1, val2;
+
+      for (uint32_t i = 0; i < colMin; i++)
+	iss >> dummy;
+      if (!(iss >> val1))
+	throw InvalidRangeException("Invalid first column");
+
+      if (col2 != col1)
+	{
+	  for (uint32_t i = 0; i < (colMax-colMin-1); i++)
+	    iss >> dummy;
+	  if (!(iss >> val2))
+	    throw InvalidRangeException("Invalid second column");
+	}
+      else
+	val2 = val2;
+
+      if (!swapped)
+	{
+	  m.a = val1;
+	  m.b = val2;
+	}
+      else
+	{
+	  m.a = val2;
+	  m.b = val1;
+	}
+
+      allData.push_back(m);
+    }
+  while (1);
+
+  sort(allData.begin(), allData.end());
+
+  double *x = new double[allData.size()];
+  double *y = new double[allData.size()];
+  
+  for (uint32_t i = 0; i < allData.size(); i++)
+    {
+      x[i] = logx ? log(allData[i].a) : allData[i].a;
+      y[i] = logy ? log(allData[i].b) : allData[i].b;
+    }
+
+  Interpolate inter = Interpolate(x, y, allData.size(), true, logx, logy);
+
+  delete[] x;
+  delete[] y;
+  
+  return inter;
+}

external/cosmotool/src/interpolate.hpp

@@ -0,0 +1,65 @@
+#ifndef __CTOOL_INTERPOLATE_HPP
+#define __CTOOL_INTERPOLATE_HPP
+
+#include "config.hpp"
+#include <inttypes.h>
+#include <gsl/gsl_interp.h>
+#include <gsl/gsl_spline.h>
+#include <vector>
+#include <utility>
+
+namespace CosmoTool
+{
+  
+  class Interpolate
+  {
+  public:
+    Interpolate() : spline(0), accel_interp(0) { }
+    Interpolate(double *xs, double *values, uint32_t N, bool autofree = false,
+		bool logx = false, bool logy = false);
+    ~Interpolate();
+
+    double compute(double x)
+      throw (InvalidRangeException);
+    double compute(double x) const
+      throw (InvalidRangeException);
+    double derivative(double x)
+      throw (InvalidRangeException);
+
+    const Interpolate& operator=(const Interpolate& a);
+
+    uint32_t getNumPoints() const;
+    void fillWithXY(double *x, double *y) const;
+    double getMaxX() const;
+    double getXi(int i) const { return spline->x[i]; }
+    double getYi(int i) const { return spline->y[i]; }
+  protected:
+    gsl_interp_accel *accel_interp;
+    gsl_spline *spline;
+    bool autoFree;
+    bool logx, logy;
+  };
+
+  typedef std::vector< std::pair<double,double> > InterpolatePairs;
+
+  Interpolate buildInterpolateFromFile(const char *fname)
+    throw (NoSuchFileException);
+  Interpolate buildInterpolateFromColumns(const char *fname, uint32_t col1, uint32_t col2, bool logx = false, bool logy = false)
+    throw (NoSuchFileException,InvalidRangeException);
+  Interpolate buildFromVector(const InterpolatePairs& v);
+
+
+  class FunctorInterpolate
+  {
+  public:
+    FunctorInterpolate(Interpolate& i) : m_i(i) {}    
+
+    double eval(double x) { return m_i.compute(x); }
+    double derivative(double x) { return m_i.derivative(x); }
+  private:
+    Interpolate& m_i;
+  };
+  
+};
+
+#endif

external/cosmotool/src/interpolate3d.hpp

@@ -0,0 +1,105 @@
+#ifndef __COSMO_INTERPOLATE3D_HPP
+#define __COSMO_INTERPOLATE3D_HPP
+
+#include "config.hpp"
+#include "field.hpp"
+#include <cmath>
+
+namespace CosmoTool
+{
+
+  template<typename IType>
+  class GridSampler
+  {
+  public:
+    typedef IType result_type;
+
+    GridSampler(IType *array_, int Nx_, int Ny_, int Nz_, int stride_)
+      : array(array_), Nx(Nx_), Ny(Ny_), Nz(Nz_), stride(stride_)
+    {
+    }
+
+    ~GridSampler()
+    {
+    }
+
+    IType& operator()(int x, int y, int z)
+      throw ()
+    {
+      while (x < 0)
+	x += Nx;
+      x %= Nx;
+      while (y < 0)
+	y += Ny;
+      y %= Ny;
+      while (z < 0)
+	z += Nz;
+      z %= Nz;
+      
+      uint32_t idx = x + Nx * (y + Ny * z);
+
+      return array[idx*stride];
+    }
+
+  private:
+    IType *array;
+    int Nx, Ny, Nz, stride;
+  };
+
+
+  // IType is the quantity to interpolate, 
+  template<typename SampledFunction, typename PosType = float>
+  class Interpolate3D
+  {
+  public:
+    typedef typename SampledFunction::result_type IType;
+
+    Interpolate3D(SampledFunction& f)
+      : sampler(f)
+    {
+    };
+
+    ~Interpolate3D()
+    {
+    };
+    
+    IType get(PosType x, PosType y, PosType z)
+      throw (InvalidArgumentException)
+    {
+      int ix = (int)std::floor(x);
+      int iy = (int)std::floor(y);
+      int iz = (int)std::floor(z);
+
+      PosType rx = x-ix;
+      PosType ry = y-iy;
+      PosType rz = z-iz;
+
+      IType v000 = sampler(ix,iy,iz);
+      IType v001 = sampler(ix,iy,iz+1);
+      IType v010 = sampler(ix,iy+1,iz);
+      IType v011 = sampler(ix,iy+1,iz+1);
+
+      IType v100 = sampler(ix+1,iy,iz);
+      IType v101 = sampler(ix+1,iy,iz+1);
+      IType v110 = sampler(ix+1,iy+1,iz);
+      IType v111 = sampler(ix+1,iy+1,iz+1);
+
+      return
+	((1-rx) * (1-ry) * (1-rz)) * v000 +
+	((1-rx) * (1-ry) *    rz)  * v001 +
+	((1-rx) *    ry  * (1-rz)) * v010 +
+	((1-rx) *    ry  *    rz)  * v011 +
+	(   rx  * (1-ry) * (1-rz)) * v100 +
+	(   rx  * (1-ry) *    rz)  * v101 +
+	(   rx  *    ry  * (1-rz)) * v110 +
+	(   rx  *    ry  *    rz)  * v111;
+    }
+
+  private:
+    SampledFunction& sampler;    
+    int Nx, Ny, Nz;
+  };
+
+};
+
+#endif

external/cosmotool/src/kdtree_leaf.hpp

@@ -0,0 +1,107 @@
+#ifndef __LEAF_KDTREE_HPP
+#define __LEAF_KDTREE_HPP
+
+#include <cmath>
+#include "config.hpp"
+#include "bqueue.hpp"
+
+namespace CosmoTool {
+
+  template<int N, typename CType = ComputePrecision> 
+  struct KDLeafDef
+  {
+    typedef CType CoordType;
+    typedef float KDLeafCoordinates[N];
+  };
+  
+  template<int N, typename ValType, typename CType = ComputePrecision>
+  struct KDLeafCell
+  {
+    bool active;
+    ValType val;
+    typename KDLeafDef<N,CType>::KDLeafCoordinates coord;
+  };
+
+  class NotEnoughCells: public Exception
+  {
+  public:
+    NotEnoughCells() : Exception() {}
+    ~NotEnoughCells() throw () {}
+  };
+
+  template<int N, typename ValType, typename CType = ComputePrecision>
+  struct KDLeafTreeNode
+  {
+    bool leaf;
+    union {
+      KDLeafCell<N,ValType,CType> *value;
+      KDLeafTreeNode<N,ValType,CType> *children[2];
+    };
+    typename KDLeafDef<N,CType>::KDLeafCoordinates minBound, maxBound;
+#ifdef __KDLEAF_TREE_NUMNODES
+    uint32_t numNodes;
+#endif
+  };
+
+  template<int N, typename ValType, typename CType = ComputePrecision>
+  class KDLeafTree
+  {
+  public:
+    typedef typename KDLeafDef<N,CType>::CoordType CoordType;
+    typedef typename KDLeafDef<N>::KDLeafCoordinates coords;
+    typedef KDLeafCell<N,ValType,CType> Cell;
+    typedef KDLeafTreeNode<N,ValType,CType> Node;
+
+    KDLeafTree(Cell *cells, uint32_t Ncells);
+    ~KDLeafTree();
+
+    Node *getRoot() { return root; }
+
+    void optimize();
+
+    Node *getAllNodes() { return nodes; }
+    uint32_t getNumNodes() const { return lastNode; }
+
+    uint32_t countActives() const;
+
+    CoordType computeDistance(const Cell *cell, const coords& x) const;
+
+#ifdef __KDLEAF_TREE_NUMNODES
+    uint32_t getNumberInNode(const Node *n) const { return n->numNodes; }
+#else
+    uint32_t getNumberInNode(const Node *n) const {
+      if (n->leaf)
+        return 1;
+      if (n == 0) 
+        return 0;
+      return getNumberInNode(n->children[0])+getNumberInNode(n->children[1]);
+    }
+#endif
+
+    double countInRange(CType sLo, CType sHi, Node *root1 = 0, Node *root2 = 0) const;
+
+  protected:
+    Node *nodes;
+    uint32_t numNodes, numCells;
+    uint32_t lastNode;
+
+    Node *root;
+    Cell **sortingHelper;
+
+    Node *buildTree(Cell **cell0,
+		    uint32_t NumCells,
+		    uint32_t depth,
+		    coords minBound,
+		    coords maxBound);
+    
+    double recursiveCountInRange(Node *na, Node *nb, CType sLo, CType sHi) const;
+  };
+
+  template<int N, typename T, typename CType>
+  uint32_t gatherActiveCells(KDLeafCell<N,T,CType> **cells, uint32_t numCells);
+
+};
+
+#include "kdtree_leaf.tcc"
+
+#endif

external/cosmotool/src/kdtree_leaf.tcc

@@ -0,0 +1,308 @@
+#include <cstring>
+#include <algorithm>
+#include <limits>
+#include <iostream>
+#include <cassert>
+
+namespace CosmoTool {
+
+  template<int N, typename ValType, typename CType>
+  class CellCompare
+  {
+  public:
+    CellCompare(int k)
+    {
+      rank = k;
+    }
+
+    bool operator()(const KDLeafCell<N,ValType,CType> *a, const KDLeafCell<N,ValType,CType> *b) const
+    {
+      return (a->coord[rank] < b->coord[rank]);
+    }
+  protected:
+    int rank;
+  };
+
+  template<int N, typename ValType, typename CType>
+  KDLeafTree<N,ValType,CType>::~KDLeafTree()
+  {
+  }
+
+  template<int N, typename ValType, typename CType>
+  KDLeafTree<N,ValType,CType>::KDLeafTree(Cell *cells, uint32_t Ncells)
+  {
+    numNodes = Ncells*3;
+    numCells = Ncells;
+    nodes = new Node[numNodes];
+
+    sortingHelper = new Cell *[Ncells];
+    for (uint32_t i = 0; i < Ncells; i++)
+	sortingHelper[i] = &cells[i];    
+
+    optimize();
+  }
+
+  template<int N, typename ValType, typename CType>
+  void KDLeafTree<N,ValType,CType>::optimize()
+  {
+    coords absoluteMin, absoluteMax;
+
+    std::cout << "Optimizing the tree..." << std::endl;
+    uint32_t activeCells = gatherActiveCells(sortingHelper, numCells);
+    std::cout << "  number of active cells = " << activeCells << std::endl;
+
+    lastNode = 0;
+    for (int i = 0; i < N; i++)
+      {
+	absoluteMin[i] = std::numeric_limits<typeof (absoluteMin[0])>::max();
+	absoluteMax[i] = -std::numeric_limits<typeof (absoluteMax[0])>::max();
+      }
+    // Find min and max corner
+    for (uint32_t i = 0; i < activeCells; i++)
+      {
+        KDLeafCell<N,ValType,CType> *cell = sortingHelper[i];
+
+        for (int k = 0; k < N; k++) {
+          if (cell->coord[k] < absoluteMin[k])
+            absoluteMin[k] = cell->coord[k];
+          if (cell->coord[k] > absoluteMax[k])
+            absoluteMax[k] = cell->coord[k];
+        }
+      }
+    
+    std::cout << "   rebuilding the tree..." << std::endl;
+    root = buildTree(sortingHelper, activeCells, 0, absoluteMin, absoluteMax);    
+    std::cout << "   done." << std::endl;
+  }
+
+  template<int N, typename ValType, typename CType>
+  uint32_t gatherActiveCells(KDLeafCell<N,ValType,CType> **cells,
+			     uint32_t Ncells)
+  {
+    uint32_t swapId = Ncells-1;
+    uint32_t i = 0;
+
+    while (!cells[swapId]->active && swapId > 0)
+      swapId--;
+
+    while (i < swapId)
+      {
+	if (!cells[i]->active)
+	  {
+	    std::swap(cells[i], cells[swapId]);
+	    while (!cells[swapId]->active && swapId > i)
+	      {
+		swapId--;
+	      }
+	  }
+	i++;
+      }
+    return swapId+1;
+  }
+
+  template<int N, typename ValType, typename CType>
+  KDLeafTreeNode<N,ValType,CType> *
+        KDLeafTree<N,ValType,CType>::buildTree(Cell **cell0,
+					       uint32_t Ncells,
+					       uint32_t depth,
+					       coords minBound,
+					       coords maxBound)
+  {
+    if (Ncells == 0)
+      return 0;
+
+    int axis = depth % N;
+    assert(lastNode != numNodes);
+    Node *node = &nodes[lastNode++];
+    uint32_t mid = Ncells/2;
+    coords tmpBound;
+    
+    // Isolate the environment
+    {
+      CellCompare<N,ValType,CType> compare(axis);
+      std::sort(cell0, cell0+Ncells, compare);
+    }
+
+    node->leaf = false;
+    memcpy(&node->minBound[0], &minBound[0], sizeof(coords));
+    memcpy(&node->maxBound[0], &maxBound[0], sizeof(coords));
+
+    if (Ncells == 1)
+      {
+	node->leaf = true;
+	node->value = *cell0;
+
+#ifdef __KDLEAF_TREE_NUMNODES
+	node->numNodes = 1;
+#endif
+	return node;
+      }
+
+    memcpy(tmpBound, maxBound, sizeof(coords));
+    tmpBound[axis] = (*(cell0+mid))->coord[axis];
+    depth++;
+    node->children[0] = buildTree(cell0, mid, depth, minBound, tmpBound);
+    
+    memcpy(tmpBound, minBound, sizeof(coords));
+    tmpBound[axis] = (*(cell0+mid))->coord[axis];
+    node->children[1] = buildTree(cell0+mid, Ncells-mid, depth, 
+				  tmpBound, maxBound);
+
+#ifdef __KDLEAF_TREE_NUMNODES
+    node->numNodes = (node->children[0] != 0) ? node->children[0]->numNodes : 0;
+    node->numNodes += (node->children[1] != 0) ? node->children[1]->numNodes : 0;
+#endif
+
+    return node;
+  }
+
+  template<int N, typename ValType, typename CType>
+  uint32_t KDLeafTree<N,ValType,CType>::countActives() const
+  {
+    uint32_t numActive = 0;
+    for (uint32_t i = 0; i < lastNode; i++)
+      {
+	if (nodes[i].value->active)
+	  numActive++;
+      }
+    return numActive;
+  }
+
+  template<int N, typename ValType, typename CType>
+  typename KDLeafDef<N,CType>::CoordType
+  KDLeafTree<N,ValType,CType>::computeDistance(const Cell *cell, const coords& x) const
+  {
+    CoordType d2 = 0;
+
+    for (int i = 0; i < N; i++)
+      {
+	CoordType delta = cell->coord[i] - x[i];
+	d2 += delta*delta;
+      }
+    return d2;
+  }
+
+ template<int N, typename ValType, typename CType>
+ double KDLeafTree<N,ValType,CType>::countInRange(CType sLo, CType sHigh, Node *root1, Node *root2) const
+ {
+   double result = recursiveCountInRange((root1 == 0) ? root : root1,
+				(root2 == 0) ? root : root2, 
+				sLo*sLo, sHigh*sHigh); 
+   return result;
+ }
+
+ 
+ template<int N, typename ValType, typename CType>
+ double KDLeafTree<N,ValType,CType>::recursiveCountInRange(Node *na, Node *nb,
+                                                         CType sLo, CType sHi) const
+ {
+   assert(nb != 0);
+   if (na == 0)
+     {
+       return 0;
+     }
+   
+   uint32_t numNa = getNumberInNode(na);
+   uint32_t numNb = getNumberInNode(nb);
+   double Cleft, Cright;
+   CType minDist, maxDist;
+
+   if (numNa == 1 && numNb == 1)
+     {
+       assert(na->leaf && nb->leaf);
+       CType ab_dist = computeDistance(na->value, nb->value->coord);
+       if (ab_dist >= sLo && ab_dist < sHi)
+	 return 1;
+       else
+	 return 0;
+     }
+   assert(numNa > 1 || numNb > 1);
+   
+   bool overlapping_a = true, overlapping_b = true;
+   for (int k = 0; k < N; k++) 
+     {
+       bool min_a_in_B = 
+	 ((na->minBound[k] >= nb->minBound[k] && 
+	   na->minBound[k] <= nb->maxBound[k]));
+       bool max_a_in_B =
+	 ((na->maxBound[k] >= nb->minBound[k] && 
+	   na->maxBound[k] <= nb->maxBound[k]));
+       bool min_b_in_A = 
+	 ((nb->minBound[k] >= na->minBound[k] && 
+	   nb->minBound[k] <= na->maxBound[k]));
+       bool max_b_in_A = 
+	 ((nb->maxBound[k] >= na->minBound[k] && 
+	   nb->maxBound[k] <= na->maxBound[k]));
+
+      if (!min_a_in_B && !max_a_in_B)
+         overlapping_a = false;
+      if (!min_b_in_A && !max_b_in_A)
+         overlapping_b = false;
+     }
+
+   if (overlapping_a || overlapping_b)
+     {
+        minDist = 0;
+        maxDist = 0;
+        for (int k = 0; k < N; k++)
+        {
+          CType delta = max(nb->maxBound[k]-na->minBound[k],na->maxBound[k]-nb->minBound[k]);
+          maxDist += delta*delta;
+        }
+     }
+   else
+     {
+       minDist = maxDist = 0;
+       for (int k = 0; k < N; k++)
+         {
+           CType delta2;
+           delta2 = max(nb->maxBound[k]-na->minBound[k],
+			na->maxBound[k]-nb->minBound[k]);
+           maxDist += delta2*delta2;
+         }
+       // mins and maxs
+       CType minmax[N][2];
+       for (int k = 0; k < N; k++)
+	 {
+	   if (na->minBound[k] < nb->minBound[k])
+	     {
+	       minmax[k][1] = na->maxBound[k];
+	       minmax[k][0] = nb->minBound[k];
+	     }
+	   else
+	     {
+	       minmax[k][1] = nb->maxBound[k];
+	       minmax[k][0] = na->minBound[k];
+	     }
+	 }
+       for (int k = 0; k < N; k++)
+	 {
+	   CType delta = max(minmax[k][0]-minmax[k][1], 0.);
+	   minDist += delta*delta;
+	 }
+     }
+   
+   if (minDist >= sHi)
+     return 0;
+   if (maxDist < sLo)
+     return 0;
+   
+   if (sLo <= minDist && maxDist < sHi)
+     return ((double)numNa)*numNb;
+   
+   if (numNa < numNb) 
+     {
+       assert(!nb->leaf);
+       Cleft = recursiveCountInRange(nb->children[0], na, sLo, sHi);
+       Cright = recursiveCountInRange(nb->children[1], na, sLo, sHi);
+     }
+   else
+     {
+       assert(!na->leaf);
+       Cleft = recursiveCountInRange(na->children[0], nb, sLo, sHi);
+       Cright = recursiveCountInRange(na->children[1], nb, sLo, sHi);
+     }
+   return Cleft+Cright;
+ }
+
+};

external/cosmotool/src/kdtree_splitters.hpp

@@ -0,0 +1,134 @@
+#ifndef __KDTREE_SPLITTERS_HPP
+#define __KDTREE_SPLITTERS_HPP
+
+#include <algorithm>
+
+namespace CosmoTool
+{
+
+  template<int N, typename ValType, typename CType = ComputePrecision> 
+  struct KD_homogeneous_cell_splitter
+  {
+    typedef typename KDDef<N,CType>::KDCoordinates coords; 
+    typedef typename KDDef<N,CType>::CoordType ctype; 
+
+
+    void check_splitting(KDCell<N,ValType,CType> **cells, uint32_t Ncells, int axis, uint32_t split_index, ctype midCoord)
+    {
+      ctype delta = std::numeric_limits<ctype>::max();
+      assert(split_index < Ncells);
+      assert(axis < N);
+      for (uint32_t i = 0; i < split_index; i++)
+	{
+	  assert(cells[i]->coord[axis] <= midCoord);
+	  delta = min(midCoord-cells[i]->coord[axis], delta);
+	}
+      for (uint32_t i = split_index+1; i < Ncells; i++)
+	{
+	  assert(cells[i]->coord[axis] > midCoord);
+	  delta = min(cells[i]->coord[axis]-midCoord, delta);
+	}
+      assert(delta >= 0);
+      assert (std::abs(cells[split_index]->coord[axis]-midCoord) <= delta);
+    } 
+    
+    void operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, coords minBound, coords maxBound)
+    {
+      if (Ncells == 1)
+	{
+	  split_index = 0;
+	  return;
+	}
+
+      ctype midCoord = 0.5*(maxBound[axis]+minBound[axis]);
+      uint32_t below = 0, above = Ncells-1;
+      ctype delta_min = std::numeric_limits<ctype>::max();
+      uint32_t idx_min = std::numeric_limits<uint32_t>::max();     
+
+      while (below < above)
+	{
+          ctype delta = cells[below]->coord[axis]-midCoord;
+	  if (delta > 0)
+            {
+              if (delta < delta_min)
+                {
+                  delta_min = delta;
+                  idx_min = above;
+                }
+	      std::swap(cells[below], cells[above--]);
+            }
+	  else
+            {
+              if (-delta < delta_min)
+                {
+                  delta_min = -delta;
+                  idx_min = below;
+                }
+	      below++;
+            }
+	}
+      // Last iteration
+      {
+	ctype delta = cells[below]->coord[axis]-midCoord;	
+	if (delta > 0)
+	  {
+	    if (delta < delta_min)
+	      {
+		delta_min = delta;
+		idx_min = above;
+	      }
+	  }
+	else
+	  {
+	    if (-delta < delta_min)
+	      {
+		delta_min = -delta;
+		idx_min = above;
+	      }
+	  }
+      }
+
+      if (idx_min != above)
+        {
+          bool cond1 = cells[idx_min]->coord[axis] > midCoord;
+          bool cond2 = cells[above]->coord[axis] > midCoord;
+          if ((cond1 && cond2) || (!cond1 && !cond2))
+            {
+              split_index = above;
+              std::swap(cells[above], cells[idx_min]);
+            }
+          else if (cond2)
+            {
+              if (above >= 1)
+		{
+		  split_index = above-1;
+		  std::swap(cells[above-1], cells[idx_min]);
+		}
+	      else
+		split_index = 0;
+	      assert(split_index >= 0);
+            }
+          else
+            {
+	      if (above+1 < Ncells)
+		{
+		  split_index = above+1;
+		  std::swap(cells[above+1], cells[idx_min]);
+		}
+	      else
+		split_index = Ncells-1;
+
+	      assert(split_index < Ncells);
+            }
+         }
+       else split_index = above;
+
+
+      //      check_splitting(cells, Ncells, axis, split_index, midCoord);
+    }
+  };
+
+
+};
+
+#endif

external/cosmotool/src/loadFlash.cpp

@@ -0,0 +1,104 @@
+/* Reads in FLASH v3 files in HDF5 format */
+
+#include <iostream>
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "load_data.hpp"
+#include "loadFlash.hpp"
+#include "h5_readFlash.hpp"
+#include "H5Cpp.h"
+
+using namespace CosmoTool;
+using namespace std;
+using namespace H5;
+
+SimuData *CosmoTool::loadFlashMulti(const char *fname, int id, int loadflags)
+{
+  SimuData *data;
+  int p, n;
+  H5File *fileID;
+  H5std_string filename;
+  //char filename[81];
+  double lbox, time, hubble, omegam, omegalambda, redshift;
+  int npart;
+
+  const double kpc2cm = 3.08568025e21;
+  const double km2cm = 1.e5;
+  const double hubble2cm = 3.240779270005e-18;
+
+  if (id != 0)
+    throw NoSuchFileException();
+
+  data = new SimuData;
+  if (data == 0) {
+    return 0;
+  }
+
+  filename = fname;
+  try {
+    H5File file (filename, H5F_ACC_RDONLY);
+
+    // simulation info
+    h5_read_flash3_header_info(&file, &time, &redshift);
+    data->time    = 1/(1+redshift);    
+
+    h5_read_runtime_parameters(&file, &lbox, &npart, &hubble, &omegam, &omegalambda);
+    data->TotalNumPart = data->NumPart = npart;
+    data->Hubble = hubble/hubble2cm;
+    data->BoxSize = lbox/kpc2cm*data->Hubble;
+    data->Omega_M = omegam;
+    data->Omega_Lambda = omegalambda;
+
+    // particle data
+    if (loadflags& NEED_POSITION) {
+      for (int i = 0; i < 3; i++) {
+        data->Pos[i] = new float[data->NumPart];
+        if (data->Pos[i] == 0) {
+	  delete data;
+	return 0;
+      }
+    } }
+
+
+    if (loadflags &NEED_VELOCITY) {
+    for (int i = 0; i < 3; i++) {
+      data->Vel[i] = new float[data->NumPart];
+      if (data->Vel[i] == 0) {
+	delete data;
+	return 0;
+      }
+    } }
+ 
+    if (loadflags & NEED_GADGET_ID) {
+    data->Id = new int[data->NumPart];
+    if (data->Id == 0) {
+      delete data;
+      return 0;
+    }
+    }
+
+    int offset = 0;
+
+    if (loadflags &(NEED_GADGET_ID|NEED_POSITION|NEED_VELOCITY))
+    h5_read_flash3_particles(&file, &npart, &npart, &offset, 
+			     data->Pos[0], data->Pos[1], data->Pos[2], 
+			     data->Vel[0], data->Vel[1], data->Vel[2],
+			     data->Id); 
+
+    for (int i = 0; i < 3; i++) {
+      for (int n = 0; n < data->NumPart; n++) {
+	if (loadflags& NEED_POSITION) data->Pos[i][n] = data->Pos[i][n] * data->Hubble / kpc2cm;
+	if (loadflags&NEED_VELOCITY) data->Vel[i][n] = data->Vel[i][n] * data->Hubble / km2cm;
+      }
+    }
+
+    file.close();
+  } catch (const FileIException& e) {
+    throw NoSuchFileException();
+  }
+
+  return data;
+}
+
+

external/cosmotool/src/loadFlash.hpp

@@ -0,0 +1,13 @@
+#ifndef __COSMO_LOAD_FLASH_HPP
+#define __COSMO_LOAD_FLASH_HPP
+
+#include "load_data.hpp"
+#include "loadSimu.hpp"
+
+namespace CosmoTool {
+  
+  SimuData *loadFlashMulti(const char *fname, int id, int flags);
+  
+};
+
+#endif

external/cosmotool/src/loadFlash_dummy.cpp

@@ -0,0 +1,9 @@
+#include "load_data.hpp"
+#include "loadFlash.hpp"
+
+using namespace CosmoTool;
+
+SimuData *CosmoTool::loadFlashMulti(const char *fname, int id, int loadflags)
+{
+  return 0;
+}

external/cosmotool/src/loadGadget.cpp

@@ -0,0 +1,349 @@
+#include <cmath>
+#include <iostream>
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "load_data.hpp"
+#include "loadGadget.hpp"
+#include "fortran.hpp"
+
+using namespace CosmoTool;
+using namespace std;
+
+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();
+  data->BoxSize = 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;
+}
+
+
+
+
+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;
+    char *out_fname = new char[k];
+    snprintf(out_fname, k, "%s.%d", fname, id);
+
+    f = new UnformattedRead(out_fname);
+    if (f == 0)
+      return 0;
+
+    delete out_fname;
+
+  } else {
+
+    f = new UnformattedRead(fname);
+    if (f == 0)
+      return 0;
+
+  }
+
+  data = new SimuData;
+  if (data == 0) {
+    delete f;
+    return 0;
+  }
+
+  long NumPart = 0, NumPartTotal = 0;
+
+  try
+    {
+      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();
+      data->time = 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();
+      data->BoxSize = h.BoxSize = f->readReal64();
+      data->Omega_M = h.Omega0 = f->readReal64();
+      data->Omega_Lambda = h.OmegaLambda = f->readReal64();
+      data->Hubble = h.HubbleParam = f->readReal64();
+      f->endCheckpoint(true);
+      
+      for(int k=0; k<6; k++)
+	{
+	  NumPart += h.npart[k];
+	  NumPartTotal += (id < 0) ? h.npart[k] : h.npartTotal[k];
+	}
+      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)
+    {
+      cerr << "Invalid format while reading header" << endl;
+      delete data;
+      delete f;
+      return 0;
+    }
+
+  if (loadflags & NEED_TYPE)
+    {
+      int p = 0;
+
+      data->type = new int[data->NumPart];
+      for (int k = 0; k < 6; k++)
+	for (int n = 0; n < h.npart[k]; n++,p++)
+	  data->type[p] = k;
+    }
+
+  if (loadflags & NEED_POSITION) {
+      
+    for (int i = 0; i < 3; i++) {
+      data->Pos[i] = new float[data->NumPart];
+      if (data->Pos[i] == 0) {
+	delete data;
+	return 0;
+      }
+    }
+    
+    try
+      {
+	f->beginCheckpoint();
+	for(int k = 0, p = 0; k < 6; k++) {
+	  for(int n = 0; n < h.npart[k]; n++) {
+	    data->Pos[0][p] = f->readReal32();
+	    data->Pos[1][p] = f->readReal32();
+	    data->Pos[2][p] = f->readReal32();
+	    p++;
+	  }
+	}
+	f->endCheckpoint();
+      }
+    catch (const InvalidUnformattedAccess& e)
+      {
+	cerr << "Invalid format while reading positions" << endl;
+	delete f;
+	delete data;
+	return 0;
+      }
+    
+  } else {
+    // Skip positions
+    f->skip(NumPart * 3 * sizeof(float) + 2*4);
+  }
+  
+  if (loadflags & NEED_VELOCITY) {
+    for (int i = 0; i < 3; i++)
+      {
+	data->Vel[i] = new float[data->NumPart];
+	if (data->Vel[i] == 0) 
+	  {
+	    delete f;
+	    delete data;
+	    return 0;
+	  }
+      }
+    
+    try
+      {
+	f->beginCheckpoint();
+	for(int k = 0, p = 0; k < 6; k++) {
+	  for(int n = 0; n < h.npart[k]; n++) {
+	    // THIS IS GADGET 1
+	    data->Vel[0][p] = f->readReal32()*velmul;
+	    data->Vel[1][p] = f->readReal32()*velmul;
+	    data->Vel[2][p] = f->readReal32()*velmul;
+	    p++;
+	  }
+	}
+	f->endCheckpoint();
+      }
+    catch (const InvalidUnformattedAccess& e)
+      {
+	cerr << "Invalid format while reading velocities" << endl;
+	delete f;
+	delete data;
+	return 0;
+      }
+    
+    // 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);
+  }
+
+  // Skip ids
+  if (loadflags & NEED_GADGET_ID) {
+    try
+      {
+	f->beginCheckpoint();
+	data->Id = new int[data->NumPart];
+	if (data->Id == 0)
+	  {
+	    delete f;
+	    delete data;
+	    return 0;
+	  }
+	
+	for(int k = 0, p = 0; k < 6; k++)
+	  {
+	    for(int n = 0; n < h.npart[k]; n++)
+	      {
+		data->Id[p] = f->readInt32();
+		p++;
+	      }
+	  }
+	f->endCheckpoint();
+      }
+    catch (const InvalidUnformattedAccess& e)
+      {
+	cerr << "Invalid formatted while reading ID" << endl;
+	delete f;
+	delete data;
+	return 0;
+      }
+  } else {
+    f->skip(2*4);
+    for (int k = 0; k < 6; k++)
+      f->skip(h.npart[k]*4);
+  }
+
+  delete f;
+
+  return data;
+}
+
+
+
+void CosmoTool::writeGadget(const char *fname, SimuData *data, int GadgetFormat)
+{
+  UnformattedWrite *f;
+  int npart[6];
+  float mass[6];
+
+  if (data->Pos[0] == 0 || data->Vel[0] == 0 || data->Id == 0)
+    return;
+
+  f = new UnformattedWrite(fname);
+  if (f == 0)
+    return;
+
+  for (int i = 0; i < 6; i++)
+    {
+      npart[i] = 0;
+      mass[i] = 0;
+    }
+  mass[1] = 1.0;
+
+  npart[1] = data->NumPart;
+  
+  f->beginCheckpoint();
+  for (int i = 0; i < 6; i++)
+    f->writeInt32(npart[i]);
+  for (int i = 0; i < 6; i++)
+    f->writeReal64(mass[i]);
+
+  f->writeReal64(data->time);
+  f->writeReal64(1/data->time-1);
+  f->writeInt32(0);
+  f->writeInt32(0);
+
+  for (int i = 0; i < 6; i++)
+    f->writeInt32(npart[i]);
+  f->writeInt32(0);
+  f->writeInt32(1);
+  f->writeReal64(data->BoxSize);
+  f->writeReal64(data->Omega_M);
+  f->writeReal64(data->Omega_Lambda);
+  f->writeReal64(data->Hubble);
+  char buf[100] = { 0, };
+  f->writeOrderedBuffer(buf, 96);
+  f->endCheckpoint();
+ 
+  f->beginCheckpoint();
+  for(int n = 0; n < data->NumPart; n++) {
+    for (int k = 0; k < 3; k++)
+      f->writeReal32(data->Pos[k][n]);
+  }
+  f->endCheckpoint();
+
+  float velmul = 1.0;
+  if (GadgetFormat == 1)
+    velmul = sqrt(data->time);
+
+  f->beginCheckpoint();
+  for(int n = 0; n < data->NumPart; n++) {
+    for (int k = 0; k < 3; k++)
+      f->writeReal32(data->Vel[k][n]/velmul);
+  }
+  f->endCheckpoint();
+
+  f->beginCheckpoint();
+  for(int n = 0; n < data->NumPart; n++)
+    {
+      f->writeInt32(data->Id[n]);
+    }
+  f->endCheckpoint();
+  delete f;
+}
+
+  

external/cosmotool/src/loadGadget.hpp

@@ -0,0 +1,18 @@
+#ifndef __COSMO_LOAD_GADGET_HPP
+#define __COSMO_LOAD_GADGET_HPP
+
+#include "load_data.hpp"
+#include "loadSimu.hpp"
+
+namespace CosmoTool {
+  
+  PurePositionData *loadGadgetPosition(const char *fname);
+
+  SimuData *loadGadgetMulti(const char *fname, int id, int flags, int GadgetFormat = 1);
+ 
+  // Only single snapshot supported
+  void writeGadget(const char *fname, SimuData *data, int GadgetFormat = 1);
+ 
+};
+
+#endif

external/cosmotool/src/loadRamses.hpp

@@ -0,0 +1,18 @@
+#ifndef _LOAD_RAMSES_HPP
+#define _LOAD_RAMSES_HPP
+
+#include "load_data.hpp"
+#include "loadSimu.hpp"
+
+namespace CosmoTool {
+
+  GadgetData *loadRamses(const char *name, bool quiet = false);
+  PurePositionData *loadRamsesPosition(const char *fname, int id, bool quiet = false, bool dp = true);
+  PhaseSpaceData *loadRamsesPhase(const char *fname, int id, bool quiet = false);
+
+  PhaseSpaceDataID *loadRamsesPhase1(const char *fname, int id, int cpuid, bool dp = true, bool quiet = false);
+
+  SimuData *loadRamsesSimu(const char *basename, int id, int cpuid, bool dp, int flags);
+};
+
+#endif

external/cosmotool/src/loadSimu.hpp

@@ -0,0 +1,48 @@
+#ifndef __COSMOTOOLBOX_HPP
+#define __COSMOTOOLBOX_HPP
+
+
+namespace CosmoTool
+{
+  static const int NEED_GADGET_ID = 1;
+  static const int NEED_POSITION = 2;
+  static const int NEED_VELOCITY = 4;
+  static const int NEED_TYPE = 8;
+
+  class SimuData 
+  {
+  public:
+    float BoxSize;
+    float time;
+    float Hubble;
+
+    float Omega_M;
+    float Omega_Lambda;
+
+    long NumPart;
+    long TotalNumPart;
+    int *Id;
+    float *Pos[3];
+    float *Vel[3];
+    int *type;
+  public:
+    SimuData() : Id(0),NumPart(0),type(0)  { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; }
+    ~SimuData() 
+    {
+      for (int j = 0; j < 3; j++)
+	{
+	  if (Pos[j])
+	    delete[] Pos[j];
+	  if (Vel[j])
+	    delete[] Vel[j];
+	}
+      if (type)
+	delete[] type;
+      if (Id)
+	delete[] Id;
+    }    
+  };
+
+};
+
+#endif

external/cosmotool/src/load_data.cpp

@@ -0,0 +1,635 @@
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "load_data.hpp"
+
+using namespace CosmoTool;
+
+//#define LARGE_CONTROL
+#define LITTLE_ENDIAN
+
+#define NEW(t,n)  ((t *)malloc(sizeof(t)*n))
+#define SKIP(f) fread(&dummy,sizeof(dummy),1,f);
+#define WRITE_DUM(f) fwrite(&dummy, sizeof(dummy),1,f);
+
+static int dummy;
+
+void CosmoTool::writeGadget(GadgetData *data, const char *fname)
+{
+  FILE *f;
+  int k, n, p;
+
+  f = fopen(fname, "w");
+  if (f == NULL) {
+    fprintf(stderr, "Cannot write gadget to file %s\n", fname);
+    return;
+  }
+
+  dummy = 256;
+  WRITE_DUM(f);
+  fwrite(&data->header, sizeof(data->header), 1, f);
+  WRITE_DUM(f);
+
+  dummy = sizeof(float)*3*data->NumPart;
+  WRITE_DUM(f);
+  for(k=0,p=0;k<5;k++) {
+    for(n=0;n<data->header.npart[k];n++) {
+      fwrite(&data->particles[p].Pos[0], sizeof(float), 3, f);
+      p++;
+    }
+  }
+  WRITE_DUM(f);
+
+  dummy = sizeof(float)*3*data->NumPart;  
+  WRITE_DUM(f);
+  for(k=0,p=0;k<6;k++) {
+    for(n=0;n<data->header.npart[k];n++) {
+      fwrite(&data->particles[p].Vel[0], sizeof(float), 3, f);
+      p++;
+    }
+  }
+  WRITE_DUM(f);
+
+  dummy = sizeof(int)*data->NumPart;
+  WRITE_DUM(f);
+  for(k=0,p=0;k<6;k++)
+    {
+      for(n=0;n<data->header.npart[k];n++)
+	{
+	  fwrite(&data->particles[p].Id, sizeof(int), 1, f);
+	  p++;
+	}
+    }
+  WRITE_DUM(f);
+
+  if(data->ntot_withmasses>0) {
+    dummy = sizeof(float)*data->NumPart;
+    WRITE_DUM(f);
+  }
+  for(k=0, p=0; k<6; k++)
+    {
+      for(n=0;n<data->header.npart[k];n++)
+	{
+	  if(data->header.mass[k]==0)
+	    fwrite(&data->particles[p].Mass, sizeof(float), 1, f);
+	  p++;
+	}
+    }
+  if(data->ntot_withmasses>0)
+    WRITE_DUM(f);      
+
+  if(data->header.npart[0]>0) {
+    dummy = data->header.npart[0]*sizeof(float);
+    WRITE_DUM(f);
+    for(n=0, p=0; n<data->header.npart[0];p++,n++) {
+      fwrite(&data->particles[p].U, sizeof(float), 1, f);
+    }
+    WRITE_DUM(f);
+    
+    WRITE_DUM(f);
+    for(n=0, p=0; n<data->header.npart[0];p++,n++) {
+      fwrite(&data->particles[p].Rho, sizeof(float), 1, f);
+    }
+    WRITE_DUM(f);
+    
+    if(data->header.flag_cooling) {
+      WRITE_DUM(f);
+      for(n=0, p=0; n<data->header.npart[0];p++,n++) {
+	fwrite(&data->particles[p].Ne, sizeof(float), 1, f);
+      }
+      WRITE_DUM(f);
+    }
+  }
+  
+  fclose(f);
+}
+
+GadgetData *CosmoTool::loadGadget(const char *fname)
+{
+  FILE *f;
+  GadgetData *data;
+  int p, k, n;
+
+  f = fopen(fname, "r");
+  if (f == NULL)
+    return NULL;
+  
+  data = NEW(GadgetData, 1);
+  SKIP(f);
+  fread(&data->header, sizeof(data->header), 1, f);
+  SKIP(f);
+  
+  for(k=0, data->ntot_withmasses=0; k<5; k++) {
+    if(data->header.mass[k]==0)
+      data->ntot_withmasses+= data->header.npart[k];
+  }
+
+  for(k=0, data->NumPart=0; k<5; k++)
+    data->NumPart+= data->header.npart[k];
+  
+  data->particles = NEW(ParticleState, data->NumPart);
+
+  SKIP(f);
+  for(k=0,p=0;k<5;k++) {
+    for(n=0;n<data->header.npart[k];n++) {
+      fread(&data->particles[p].Pos[0], sizeof(float), 3, f);
+      p++;
+    }
+  }
+  SKIP(f);
+  
+  SKIP(f);
+  for(k=0,p=0;k<6;k++) {
+    for(n=0;n<data->header.npart[k];n++) {
+      fread(&data->particles[p].Vel[0], sizeof(float), 3, f);
+      p++;
+    }
+  }
+  SKIP(f);
+  
+
+  SKIP(f);
+  for(k=0,p=0;k<6;k++)
+    {
+      for(n=0;n<data->header.npart[k];n++)
+	{
+	  fread(&data->particles[p].Id, sizeof(int), 1, f);
+	  p++;
+	}
+    }
+  SKIP(f);
+
+  if(data->ntot_withmasses>0)
+    SKIP(f);
+  for(k=0, p=0; k<6; k++)
+    {
+      for(n=0;n<data->header.npart[k];n++)
+	{
+	  data->particles[p].Type=k;
+
+	  if(data->header.mass[k]==0)
+	    fread(&data->particles[p].Mass, sizeof(float), 1, f);
+	  else
+	    data->particles[p].Mass= data->header.mass[k];
+	  p++;
+	}
+    }
+  if(data->ntot_withmasses>0)
+    SKIP(f);      
+
+  if(data->header.npart[0]>0)
+    {
+      SKIP(f);
+      for(n=0, p=0; n<data->header.npart[0];p++,n++) {
+	  fread(&data->particles[p].U, sizeof(float), 1, f);
+	}
+      SKIP(f);
+
+      SKIP(f);
+      for(n=0, p=0; n<data->header.npart[0];p++,n++) {
+	fread(&data->particles[p].Rho, sizeof(float), 1, f);
+      }
+      SKIP(f);
+
+      if(data->header.flag_cooling)
+	{
+	  SKIP(f);
+	  for(n=0, p=0; n<data->header.npart[0];p++,n++)
+	    {
+	      fread(&data->particles[p].Ne, sizeof(float), 1, f);
+	    }
+	  SKIP(f);
+	}
+      else
+	for(n=0, p=0; n<data->header.npart[0];p++,n++)
+	  {
+	    data->particles[p].Ne= 1.0;
+	  }
+    }
+
+
+  fclose(f);
+
+  return data;
+}
+
+void CosmoTool::freeGadget(GadgetData *data)
+{
+  free(data->particles);
+  free(data);
+}
+
+void CosmoTool::writePersoSet(ParticleSet *set, const char *fname)
+{
+  FILE *f;
+  int i;
+
+  f = fopen(fname, "w");
+  if (f == NULL) {
+    perror("writePersoSet");
+    return;
+  }
+
+  fwrite(&set->header, sizeof(set->header), 1, f);
+  fwrite(set->Npart, sizeof(set->Npart[0]), set->header.Ntypes, f);
+  
+  for (i=0;i<set->header.Ntypes;i++)
+    fwrite(set->particles[i], sizeof(ParticleState), set->Npart[i], f);
+
+  fclose(f);
+}
+
+ParticleSet *CosmoTool::loadPersoSet(const char *fname)
+{
+  ParticleSet *set;
+  FILE *f;
+  int i;
+
+  f = fopen(fname, "r");
+  if (f == NULL) {
+    perror("loadPersoSet");
+    return NULL;
+  }
+
+  set = NEW(ParticleSet, 1);
+  fread(&set->header, sizeof(set->header), 1, f);
+
+  set->Npart = NEW(int, set->header.Ntypes);
+  fread(set->Npart, sizeof(set->Npart[0]), set->header.Ntypes, f);;
+  
+  set->particles = NEW(ParticleState *, set->header.Ntypes);
+  for (i=0;i<set->header.Ntypes;i++) {
+    set->particles[i] = NEW(ParticleState, set->Npart[i]);
+    fread(set->particles[i], sizeof(ParticleState), set->Npart[i], f);
+  }
+
+  fclose(f);
+
+  return set;
+}
+
+void CosmoTool::freePersoSet(ParticleSet *set)
+{
+  int i;
+
+  for (i=0;i<set->header.Ntypes;i++) {
+    free(set->particles[i]);
+  }
+  if (set->Npart != NULL) {
+    free(set->particles);
+    free(set->Npart);
+  }
+}
+
+#ifdef WANT_MAIN
+int main(int argc, char **argv) {
+  GadgetData *data;
+  FILE *plot;
+  int i;
+  double bl;
+  int N;
+  double rms;
+
+  if (argc < 3) {
+      fprintf(stderr, "Usage: %s [GADGET DATA FILE] [BOXSIZE] [N PARTIC]\n", argv[0]);
+      return -1;
+  }
+
+  plot = fopen("plot", "w");
+  
+  bl = atof(argv[2]);
+  data = loadGadget(argv[1]);
+
+  printf("Redshift: %lg\n", data->header.redshift);
+  rms = 0;
+  N = atoi(argv[3]);
+  for (i=0;i<data->NumPart;i++) {
+    if (i == data->header.npart[0])
+	    fprintf(plot,"\n\n");
+
+    fprintf(plot, "%f %f %f\n", data->particles[i].Pos[0], data->particles[i].Pos[1], data->particles[i].Pos[2]);
+
+
+    /* Compute the RMS */
+    {
+      /* First find the nearest grid node. */
+      int k;
+      int x;
+      double dx;
+
+      for (k=0;k<3;k++) {
+	x = data->particles[i].Pos[k] / bl * N;
+	dx = data->particles[i].Pos[k]-x*bl/N;
+	rms += dx*dx;
+      }
+    }
+  }
+
+  printf("delta rms = %e\n", sqrt(rms/data->NumPart));
+  freeGadget(data);
+
+  fclose(plot);
+
+  return 0;
+}
+#endif
+
+#define LEN0 200.0
+
+GadgetData *CosmoTool::loadSimulationData(const char *fname)
+{
+  GadgetData *gd = NEW(GadgetData, 1);
+  FILE *f;
+  int lineNo;
+  char line[1024];
+  int i;
+  int j;
+  
+  gd->header.BoxSize = LEN0;
+
+  f = fopen(fname, "r");
+  lineNo = 0;
+  while (!feof(f))
+    {
+      fgets(line, sizeof(line), f);
+      lineNo++;
+    }
+  lineNo--;
+  rewind(f);
+
+  gd->NumPart = lineNo;
+  gd->particles = NEW(ParticleState, lineNo);
+
+  i = 0;
+  while (!feof(f))
+    {
+      fgets(line, sizeof(line), f);
+      int r = sscanf(line, "%*d %*d %f %f %f %f %f %f %f %f %f %*f %*f %*f %f %f %f",
+		     &gd->particles[i].Pos[0], &gd->particles[i].Pos[1], &gd->particles[i].Pos[2],
+		     &gd->particles[i].Init[0], &gd->particles[i].Init[1], &gd->particles[i].Init[2],
+		     &gd->particles[i].Vel[0], &gd->particles[i].Vel[1], &gd->particles[i].Vel[2],
+		     &gd->particles[i].VelInit[0], &gd->particles[i].VelInit[1], &gd->particles[i].VelInit[2]
+		     );
+      if (r != 12)
+	{
+	  printf("line %d: '%s'\n", i, line);
+	  printf("returned r=%d\n", r);
+	  abort();
+	}
+      assert(r == 12);
+      for (j = 0; j < 3; j++)
+	{
+	  gd->particles[i].Vel[j] *= 100.0 * LEN0 / (0.9641010);
+	  gd->particles[i].VelInit[j] *= 100.0 * 1/71. * LEN0 / (0.9641010);
+	  gd->particles[i].Pos[j] *= LEN0;
+	  gd->particles[i].Init[j] *= LEN0;
+	}
+
+      gd->particles[i].Type = 0;
+      gd->particles[i].Mass = 1.0;
+      gd->particles[i].Id = i;
+
+      i++;
+    }
+  fclose(f);
+
+  return gd;
+}
+
+#ifndef LITTLE_ENDIAN 
+#define read_buf(b, n) \
+{ \
+  int k; \
+  control_size -= n; \
+  for (k = (n-1); k >= 0; k--) \
+    fread(&b[k], 1, 1, infile); \
+}
+#else
+#define read_buf(b, n) \
+{ \
+  int k; \
+  control_size -= n; \
+  for (k = 0; k < n; k++) \
+    fread(&b[k], 1, 1, infile); \
+}
+#endif
+
+#define read_int(i) \
+{ \
+  char *o = (char*)&(i); \
+  read_buf(o, 4); \
+}
+
+#define read_real(f) \
+{ \
+  char *o = (char*)&(f); \
+  read_buf(o, 4); \
+}
+
+#define read_characters(c, n) { \
+  int k; \
+  control_size -= n; \
+  fread(c, 1, n, outfile); \
+}
+
+#define push_dummy_control(id) \
+{ int control_size = 0;
+
+#define pop_dummy_control() }
+
+#if defined(LARGE_CONTROL) && defined(LITTLE_ENDIAN)
+#define push_control(id) \
+{ \
+  int control_size = 0; \
+  int control_size2 = 0; \
+  char *intbuf = (char*)&control_size; \
+  fread(&control_size, 8, 1, infile);
+
+#define pop_control(id) \
+  fread(&control_size2, 8, 1, infile); \
+  assert(control_size == 0); \
+}
+#elif !defined(LARGE_CONTROL) && defined(LITTLE_ENDIAN)
+#define push_control(id) \
+{ \
+  int control_size = 0; \
+  int control_size2 = 0; \
+  char *intbuf = (char*)&control_size; \
+  fread(&control_size, 4, 1, infile);
+
+#define pop_control(id) \
+  fread(&control_size2, 4, 1, infile); \
+  assert(control_size == 0); \
+}
+
+#elif defined(LARGE_CONTROL) && !defined(LITTLE_ENDIAN)
+#define push_control(id) \
+{ \
+  int control_size = 0; \
+  int control_size2 = 0; \
+  char *intbuf = (char*)&control_size; \
+  fread(&control_size, 8, 1, infile);
+
+#define pop_control(id) \
+  fread(&control_size2, 8, 1, infile); \
+  assert(control_size == 0); \
+}
+
+#elif !defined(LARGE_CONTROL) && !defined(LITTLE_ENDIAN)
+#define push_control(id) \
+{ \
+  int control_size = 0; \
+  int control_size2 = 0; \
+  char *intbuf = (char*)&control_size; \
+  fread(&control_size, 4, 1, infile);
+
+#define pop_control(id) \
+  fread(&control_size2, 4, 1, infile); \
+  assert(control_size == 0); \
+}
+
+#endif
+
+GadgetData *CosmoTool::loadHydra(const char *fname)
+{
+  GadgetData *gd = NEW(GadgetData, 1);
+  FILE *f;
+  int version0, version1, version2;
+  int irun, nobj, ngas, ndark, intl, nlmx, perr;
+  float dtnorm, sft0, sftmin, sftmax;
+  int pad3;
+  float h100, box100, zmet0;
+  int lcool;
+  float rmnorm0;
+  int pad4, pad5;
+  float tstart, omega0, xlambda0, h0t0, rcen, rmax2;
+  float rmbary;
+  int j;
+  float atime;
+
+  f = fopen(fname, "r");
+#define infile f
+  push_control(0);
+  read_int(version0);
+  read_int(version1);
+  read_int(version2);
+  pop_control(0);
+
+  if (version0 != 4)
+    {
+      fclose(f);
+      return NULL;
+    }
+  push_control(1);
+  
+  for (j = 0; j < 200; j++)
+    {
+      int mydummy;
+      read_int(mydummy);
+    }
+  for (j = 0; j < 5; j++)
+  {
+	  float mydummy;
+	  read_real(mydummy);
+  }
+  read_real(atime);
+  gd->header.time = atime;
+  gd->header.redshift = 1/atime - 1;
+
+  for (j = 6; j < 100; j++)
+    {
+      int mydummy;
+      read_int(mydummy);
+    }
+  read_int(irun);
+  read_int(nobj);
+  read_int(ngas);
+  read_int(ndark);
+  read_int(intl);
+  read_int(nlmx);
+  read_int(perr);
+  read_real(dtnorm);
+  read_real(sft0);
+  read_real(sftmin);
+  read_real(sftmax);
+  read_int(pad3);
+  read_real(h100);
+  read_real(box100);
+  read_real(zmet0);
+  read_int(lcool);
+  read_real(rmbary);
+  read_real(rmnorm0);
+  read_int(pad4);
+  read_int(pad5);
+  read_real(tstart);
+  read_real(omega0);
+  read_real(xlambda0);
+  read_real(h0t0);
+  read_real(rcen);
+  read_real(rmax2);
+  for (j = 0; j < 74; j++)
+    {
+      int mydummy;
+      read_int(mydummy);
+    }
+  pop_control(1);
+
+  gd->header.npart[1] = ndark;
+  gd->header.npart[0] = ngas;
+  gd->header.num_files = 1;
+  gd->header.flag_cooling = lcool;
+  gd->header.BoxSize = box100 * 1000;
+  gd->header.HubbleParam = h100;
+  gd->header.Omega0 = omega0;
+  gd->header.OmegaLambda = xlambda0;
+
+  push_control(2);
+  for (j = 0; j < nobj; j++)
+    {
+      int mydummy;
+      read_int(mydummy);
+    }
+  pop_control(2);
+
+  gd->NumPart = nobj;
+  gd->ntot_withmasses = nobj;
+  gd->particles = NEW(ParticleState, nobj);
+
+  push_control(3);
+  for (j = 0; j < nobj; j++)
+    {
+      float rm;
+      gd->particles[j].Id = j;
+      read_real(gd->particles[j].Mass);
+    }
+  pop_control(3);
+
+  push_control(4);
+  for (j = 0; j < nobj; j++)
+    {
+      int k;
+      for (k = 0; k < 3; k++)
+	{
+	  read_real(gd->particles[j].Pos[k]);
+	  gd->particles[j].Pos[k] *= gd->header.BoxSize;
+	}
+    }
+  pop_control(4);
+
+  push_control(5);
+  for (j = 0; j < nobj; j++)
+    {
+      int k;
+      for (k = 0; k < 3; k++)
+	{
+	  read_real(gd->particles[j].Vel[k]);
+	  gd->particles[j].Vel[k] *= 100.0 * box100  / h0t0 * atime;
+	}
+    }
+  pop_control(5);
+
+  fclose(f);
+#undef infile
+
+  return gd;
+}

external/cosmotool/src/load_data.hpp

@@ -0,0 +1,97 @@
+#ifndef _LOAD_GADGET_DATA_HPP
+#define _LOAD_GADGET_DATA_HPP
+
+#include "config.hpp"
+
+namespace CosmoTool {
+
+  struct GadgetHeader
+  {
+    int      npart[6];
+    double   mass[6];
+    double   time;
+    double   redshift;
+    int      flag_sfr;
+    int      flag_feedback;
+    int      npartTotal[6];
+    int      flag_cooling;
+    int      num_files;
+    double   BoxSize;
+    double   Omega0;
+    double   OmegaLambda;
+    double   HubbleParam; 
+    char     fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8];  /* fills to 256 Bytes */
+  };
+
+  struct ParticleState
+  {
+    float  Pos[3];
+    float  Init[3];
+    float  Vel[3];
+    float  VelInit[3];
+    float  Mass;
+    int    Type;
+    
+    float  Rho, U, Temp, Ne;
+    int Id;
+  };
+  
+  struct GadgetData {
+    GadgetHeader header;
+    ParticleState *particles;
+    int NumPart;
+    int ntot_withmasses;
+  };
+  
+  struct ParticleSetHeader {
+    int Ntypes;
+    float BoxSize;
+    float RedShift;
+    char header[256 - 4 - 2*4];
+  };
+
+  struct ParticleSet {
+    ParticleSetHeader header;
+    // Particle description
+    int *Npart;
+    ParticleState **particles;
+  };
+
+  struct PurePositionData {
+    unsigned int NumPart;
+    double BoxSize;
+    double hubble;
+    FCoordinates *pos;
+  };
+
+  struct PhaseSpaceData {
+    unsigned int NumPart;
+    double hubble;
+    double BoxSize;
+    FCoordinates *pos;
+    FCoordinates *vel;
+  };
+
+  struct PhaseSpaceDataID {
+    unsigned int NumPart;
+    double hubble;
+    double BoxSize;
+    FCoordinates *pos;
+    FCoordinates *vel;
+    int *ID;
+  };
+
+  void writeGadget(GadgetData *data, const char *fname);
+  GadgetData *loadGadget(const char *fname);
+  void freeGadget(GadgetData *data);
+
+  GadgetData *loadSimulationData(const char *fname);
+  GadgetData *loadHydra(const char *fname);
+
+  void writePersoSet(ParticleSet *set, const char *fname);
+  ParticleSet *loadPersoSet(const char *fname);
+  void freePersoSet(ParticleSet *set);
+
+};
+
+#endif

external/cosmotool/src/mach.hpp

@@ -0,0 +1,20 @@
+#ifndef __COSMO_MACHINE_TEST_HPP
+#define __COSMO_MACHINE_TEST_HPP
+
+#include <iostream>
+
+template<typename T>
+T mach_epsilon()
+{
+  T eps = (T)1;
+
+  do
+    {
+      eps /= 2;
+    }
+  while ((T)(1 + (eps/2)) != (T)1);
+
+  return eps;
+}
+
+#endif

external/cosmotool/src/miniargs.cpp

@@ -0,0 +1,55 @@
+#include <cstdlib>
+#include <cstdio>
+#include <cstring>
+#include "miniargs.hpp"
+#include <iostream>
+
+using namespace CosmoTool;
+using namespace std;
+
+int CosmoTool::parseMiniArgs(int argc, char **argv, MiniArgDesc *desc)
+{
+  int numMini;
+  for (numMini = 0; desc[numMini].name != 0; numMini++);
+
+  if ((argc-1) != numMini)
+    {
+      cerr << "Usage: ";
+      for (int i = 0; i < numMini; i++)
+	{
+	  cerr << '[' << desc[i].name << "] ";
+	}
+      cerr << endl;
+      return 0;
+    }
+
+  for (int i = 0; i < numMini; i++)
+    {
+      switch (desc[i].argType)
+	{
+	case MINIARG_STRING:
+	  *((char **)desc[i].data) = strdup(argv[i+1]);
+	  break;
+	case MINIARG_INT:
+	  *((int *)desc[i].data) = strtol(argv[i+1], NULL, 0);
+	  break;
+	case MINIARG_DOUBLE:
+	  *((double *)desc[i].data) = strtod(argv[i+1], NULL);
+	  break;
+        case MINIARG_FLOAT:
+	  *((float *)desc[i].data) = strtod(argv[i+1], NULL);
+          break;
+	case MINIARG_DOUBLE_3D_VECTOR:
+	  {
+	    double *d_array = (double *)(desc[i].data);
+
+	    if (sscanf(argv[i+1], "(%lg,%lg,%lg)", &d_array[0], &d_array[1], &d_array[2]) != 3)
+	      return 0;
+	    break;
+	  }
+	}
+    }
+
+  return 1;
+}
+

external/cosmotool/src/miniargs.hpp

@@ -0,0 +1,26 @@
+#ifndef _MAK_MINIARGS_HPP
+#define _MAK_MINIARGS_HPP
+
+namespace CosmoTool
+{
+  typedef enum 
+    {
+      MINIARG_NULL,
+      MINIARG_STRING,
+      MINIARG_INT,
+      MINIARG_DOUBLE,
+      MINIARG_FLOAT,
+      MINIARG_DOUBLE_3D_VECTOR
+    } MiniArgType;
+  
+  typedef struct
+  {
+    const char *name;
+    void *data;
+    MiniArgType argType;
+  } MiniArgDesc;
+
+  int parseMiniArgs(int argc, char **argv, MiniArgDesc *desc);
+};
+
+#endif

external/cosmotool/src/mykdtree.hpp

@@ -0,0 +1,183 @@
+#ifndef __HV_KDTREE_HPP
+#define __HV_KDTREE_HPP
+
+#include <cmath>
+#include "config.hpp"
+#include "bqueue.hpp"
+
+namespace CosmoTool {
+
+  template<int N, typename CType = ComputePrecision> 
+  struct KDDef
+  {
+    typedef CType CoordType;
+    typedef float KDCoordinates[N];
+  };
+  
+  template<int N, typename ValType, typename CType = ComputePrecision>
+  struct KDCell
+  {
+    bool active;
+    ValType val;
+    typename KDDef<N,CType>::KDCoordinates coord;
+  };
+
+  class NotEnoughCells: public Exception
+  {
+  public:
+    NotEnoughCells() : Exception() {}
+    ~NotEnoughCells() throw () {}
+  };
+
+  class InvalidOnDiskKDTree : public Exception
+  {
+  public:
+    InvalidOnDiskKDTree() : Exception() {}
+    ~InvalidOnDiskKDTree() throw () {}
+  };
+
+  template<int N, typename ValType, typename CType = ComputePrecision>
+  struct KDTreeNode
+  {
+    KDCell<N,ValType,CType> *value;
+    KDTreeNode<N,ValType,CType> *children[2];
+    typename KDDef<N,CType>::KDCoordinates minBound, maxBound;
+#ifdef __KD_TREE_NUMNODES
+    uint32_t numNodes;
+#endif
+  };
+
+  template<int N, typename ValType, typename CType = ComputePrecision>
+  class RecursionInfoCells
+  {
+  public:
+
+    typename KDDef<N,CType>::KDCoordinates x;
+    typename KDDef<N,CType>::CoordType r, r2;
+    KDCell<N, ValType,CType> **cells;
+    typename KDDef<N,CType>::CoordType *distances;
+    uint32_t currentRank;
+    uint32_t numCells;
+  };
+  
+
+  template<int N, typename ValType, typename CType = ComputePrecision>
+  class RecursionMultipleInfo
+  {
+  public:
+    const typename KDDef<N,CType>::KDCoordinates& x;
+    BoundedQueue< KDCell<N,ValType,CType> *, typename KDDef<N,CType>::CoordType> queue;
+    int traversed;
+
+    RecursionMultipleInfo(const typename KDDef<N,CType>::KDCoordinates& rx,
+			  KDCell<N,ValType,CType> **cells,
+			  uint32_t numCells)
+      : x(rx), queue(cells, numCells, INFINITY),traversed(0)
+    {      
+    }
+  };
+
+  template<int N, typename ValType, typename CType = ComputePrecision> 
+  struct KD_default_cell_splitter
+  {
+    void operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound);
+  };
+
+  template<int N, typename ValType, typename CType = ComputePrecision, typename CellSplitter = KD_default_cell_splitter<N,ValType,CType> >
+  class KDTree
+  {
+  public:
+    typedef typename KDDef<N,CType>::CoordType CoordType;
+    typedef typename KDDef<N>::KDCoordinates coords;
+    typedef KDCell<N,ValType,CType> Cell;
+    typedef KDTreeNode<N,ValType,CType> Node;
+    
+    CellSplitter splitter;
+
+    KDTree(Cell *cells, uint32_t Ncells);
+    ~KDTree();
+
+    uint32_t getIntersection(const coords& x, CoordType r, 
+			     Cell **cells,
+			     uint32_t numCells)
+      throw (NotEnoughCells);
+    uint32_t getIntersection(const coords& x, CoordType r, 
+			     Cell **cells,
+			     CoordType *distances,
+			     uint32_t numCells)
+      throw (NotEnoughCells);
+    uint32_t countCells(const coords& x, CoordType r);
+
+    Cell *getNearestNeighbour(const coords& x);
+
+    void getNearestNeighbours(const coords& x, uint32_t NumCells,
+			      Cell **cells);
+    void getNearestNeighbours(const coords& x, uint32_t NumCells,
+			      Cell **cells,
+			      CoordType *distances);
+
+    Node *getRoot() { return root; }
+
+    void optimize();
+
+    Node *getAllNodes() { return nodes; }
+    uint32_t getNumNodes() const { return lastNode; }
+
+    uint32_t countActives() const;
+
+#ifdef __KD_TREE_NUMNODES
+    uint32_t getNumberInNode(const Node *n) const { return n->numNodes; }
+#else
+    uint32_t getNumberInNode(const Node *n) const {
+      if (n == 0) 
+        return 0;
+      return 1+getNumberInNode(n->children[0])+getNumberInNode(n->children[1]);
+    }
+#endif
+
+#ifdef __KD_TREE_SAVE_ON_DISK
+    KDTree(std::istream& i, Cell *cells, uint32_t Ncells)
+      throw (InvalidOnDiskKDTree);
+
+    void saveTree(std::ostream& o) const;
+#endif
+  protected:
+    Node *nodes;
+    uint32_t numNodes;
+    uint32_t lastNode;
+
+    Node *root;
+    Cell **sortingHelper;
+    Cell *base_cell;
+
+    Node *buildTree(Cell **cell0,
+		    uint32_t NumCells,
+		    uint32_t depth,
+		    coords minBound,
+		    coords maxBound);
+    
+    template<bool justCount>
+    void recursiveIntersectionCells(RecursionInfoCells<N,ValType, CType>& info,
+				    Node *node,
+				    int level)
+      throw (NotEnoughCells);
+
+    CoordType computeDistance(const Cell *cell, const coords& x) const;
+    void recursiveNearest(Node *node,
+			  int level,
+			  const coords& x,
+			  CoordType& R2,
+			  Cell*& cell);
+    void recursiveMultipleNearest(RecursionMultipleInfo<N,ValType,CType>& info, Node *node,
+				  int level);    
+
+  };
+
+  template<int N, typename T, typename CType>
+  uint32_t gatherActiveCells(KDCell<N,T,CType> **cells, uint32_t numCells);
+
+};
+
+#include "mykdtree.tcc"
+
+#endif

external/cosmotool/src/mykdtree.tcc

@@ -0,0 +1,599 @@
+#include <cstring>
+#include <algorithm>
+#include <limits>
+#include <iostream>
+#include <cassert>
+
+namespace CosmoTool {
+
+  template<int N, typename ValType, typename CType>
+  class CellCompare
+  {
+  public:
+    CellCompare(int k)
+    {
+      rank = k;
+    }
+
+    bool operator()(const KDCell<N,ValType,CType> *a, const KDCell<N,ValType,CType> *b) const
+    {
+      return (a->coord[rank] < b->coord[rank]);
+    }
+  protected:
+    int rank;
+  };
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  KDTree<N,ValType,CType,CellSplitter>::~KDTree()
+  {
+  }
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  KDTree<N,ValType,CType,CellSplitter>::KDTree(Cell *cells, uint32_t Ncells)
+  {
+    
+    base_cell = cells;
+    numNodes = Ncells;
+    nodes = new Node[numNodes];
+
+    sortingHelper = new Cell *[Ncells];
+    for (uint32_t i = 0; i < Ncells; i++)
+	sortingHelper[i] = &cells[i];    
+
+    optimize();
+  }
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  void KDTree<N,ValType,CType,CellSplitter>::optimize()
+  {
+    coords absoluteMin, absoluteMax;
+
+    std::cout << "Optimizing the tree..." << std::endl;
+    uint32_t activeCells = gatherActiveCells(sortingHelper, numNodes);
+    std::cout << "  number of active cells = " << activeCells << std::endl;
+
+    lastNode = 0;
+    for (int i = 0; i < N; i++)
+      {
+	absoluteMin[i] = std::numeric_limits<typeof (absoluteMin[0])>::max();
+	absoluteMax[i] = -std::numeric_limits<typeof (absoluteMax[0])>::max();
+      }
+    // Find min and max corner
+    for (uint32_t i = 0; i < activeCells; i++)
+      {
+        KDCell<N,ValType,CType> *cell = sortingHelper[i];
+
+        for (int k = 0; k < N; k++) {
+          if (cell->coord[k] < absoluteMin[k])
+            absoluteMin[k] = cell->coord[k];
+          if (cell->coord[k] > absoluteMax[k])
+            absoluteMax[k] = cell->coord[k];
+        }
+      }
+    
+    std::cout << "   rebuilding the tree..." << std::endl;
+    root = buildTree(sortingHelper, activeCells, 0, absoluteMin, absoluteMax);    
+    std::cout << "   done." << std::endl;
+  }
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r, 
+						    KDTree<N,ValType,CType,CellSplitter>::Cell **cells,
+						    uint32_t numCells)
+    throw (NotEnoughCells)
+  {
+    RecursionInfoCells<N,ValType,CType> info;
+
+    memcpy(info.x, x, sizeof(x));
+    info.r = r;
+    info.r2 = r*r;
+    info.cells = cells;
+    info.currentRank = 0;
+    info.numCells = numCells;
+    info.distances = 0;
+
+    recursiveIntersectionCells<false>(info, root, 0);
+    return info.currentRank;
+  }
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r, 
+						    Cell **cells,
+						    CoordType *distances,
+						    uint32_t numCells)
+    throw (NotEnoughCells)
+  {
+    RecursionInfoCells<N,ValType> info;
+
+    memcpy(info.x, x, sizeof(x));
+    info.r = r;
+    info.r2 = r*r;
+    info.cells = cells;
+    info.currentRank = 0;
+    info.numCells = numCells;
+    info.distances = distances;
+
+    recursiveIntersectionCells<false>(info, root, 0);
+    return info.currentRank;
+  }
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  uint32_t KDTree<N,ValType,CType,CellSplitter>::countCells(const coords& x, CoordType r)
+  {
+    RecursionInfoCells<N,ValType> info;
+
+    memcpy(info.x, x, sizeof(x));
+    info.r = r;
+    info.r2 = r*r;
+    info.cells = 0;
+    info.currentRank = 0;
+    info.numCells = 0;
+    info.distances = 0;
+
+    recursiveIntersectionCells<true>(info, root, 0);
+    return info.currentRank;
+  }
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  template<bool justCount>
+  void KDTree<N,ValType,CType,CellSplitter>::recursiveIntersectionCells(RecursionInfoCells<N,ValType,CType>& info, 
+							   Node *node,
+							   int level)
+    throw (NotEnoughCells)
+  {
+    int axis = level % N;
+    CoordType d2 = 0;
+
+    if (node->value->active) 
+      {
+	for (int j = 0; j < 3; j++)
+	  {
+	    CoordType delta = info.x[j]-node->value->coord[j];
+	    d2 += delta*delta;
+	  }
+	if (d2 < info.r2)
+	  {
+	    if (!justCount)
+	      {
+		if (info.currentRank == info.numCells)
+		  throw NotEnoughCells();
+		info.cells[info.currentRank] = node->value;
+		if (info.distances)
+		  info.distances[info.currentRank] = d2;
+	      }
+	    info.currentRank++;
+	  }
+      }
+
+    // The hypersphere intersects the left child node
+    if (((info.x[axis]+info.r) > node->minBound[axis]) &&
+	((info.x[axis]-info.r) < node->value->coord[axis]))
+      {
+	if (node->children[0] != 0)
+	  recursiveIntersectionCells<justCount>(info, node->children[0],
+				     level+1);
+      }
+    if (((info.x[axis]+info.r) > node->value->coord[axis]) &&
+	((info.x[axis]-info.r) < node->maxBound[axis]))
+      {
+	if (node->children[1] != 0)
+	  recursiveIntersectionCells<justCount>(info, node->children[1],
+					       level+1);
+      }
+  }
+
+  template<int N, typename ValType, typename CType>
+  uint32_t gatherActiveCells(KDCell<N,ValType,CType> **cells,
+			     uint32_t Ncells)
+  {
+    uint32_t swapId = Ncells-1;
+    uint32_t i = 0;
+
+    while (!cells[swapId]->active && swapId > 0)
+      swapId--;
+
+    while (i < swapId)
+      {
+	if (!cells[i]->active)
+	  {
+	    std::swap(cells[i], cells[swapId]);
+	    while (!cells[swapId]->active && swapId > i)
+	      {
+		swapId--;
+	      }
+	  }
+	i++;
+      }
+    return swapId+1;
+  }
+
+  template<int N, typename ValType, typename CType>
+  void  KD_default_cell_splitter<N,ValType,CType>::operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound)
+  {
+    CellCompare<N,ValType,CType> compare(axis);
+    std::sort(cells, cells+Ncells, compare);
+    split_index = Ncells/2;
+  }
+
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  KDTreeNode<N,ValType,CType> *KDTree<N,ValType,CType,CellSplitter>::buildTree(Cell **cell0,
+								  uint32_t Ncells,
+								  uint32_t depth,
+								  coords minBound,
+								  coords maxBound)
+  {
+    if (Ncells == 0)
+      return 0;
+
+    int axis = depth % N;
+    Node *node = &nodes[lastNode++];
+    uint32_t mid;
+    coords tmpBound;
+    
+    // Isolate the environment
+    splitter(cell0, Ncells, mid, axis, minBound, maxBound);
+
+    node->value = *(cell0+mid);
+    memcpy(&node->minBound[0], &minBound[0], sizeof(coords));
+    memcpy(&node->maxBound[0], &maxBound[0], sizeof(coords));
+
+    memcpy(tmpBound, maxBound, sizeof(coords));
+    tmpBound[axis] = node->value->coord[axis];
+
+    depth++;
+    node->children[0] = buildTree(cell0, mid, depth, minBound, tmpBound);
+    
+    memcpy(tmpBound, minBound, sizeof(coords));
+    tmpBound[axis] = node->value->coord[axis];
+    node->children[1] = buildTree(cell0+mid+1, Ncells-mid-1, depth, 
+				  tmpBound, maxBound);
+
+#ifdef __KD_TREE_NUMNODES
+    node->numNodes = (node->children[0] != 0) ? node->children[0]->numNodes : 0;
+    node->numNodes += (node->children[1] != 0) ? node->children[1]->numNodes : 0;
+    node->numNodes++;
+#endif
+
+    return node;
+  }
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  uint32_t KDTree<N,ValType,CType,CellSplitter>::countActives() const
+  {
+    uint32_t numActive = 0;
+    for (uint32_t i = 0; i < lastNode; i++)
+      {
+	if (nodes[i].value->active)
+	  numActive++;
+      }
+    return numActive;
+  }
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  typename KDDef<N,CType>::CoordType
+  KDTree<N,ValType,CType,CellSplitter>::computeDistance(const Cell *cell, const coords& x) const
+  {
+    CoordType d2 = 0;
+
+    for (int i = 0; i < N; i++)
+      {
+	CoordType delta = cell->coord[i] - x[i];
+	d2 += delta*delta;
+      }
+    return d2;
+  }
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  void
+  KDTree<N,ValType,CType,CellSplitter>::recursiveNearest(
+				       Node *node,
+				       int level,
+				       const coords& x,
+				       CoordType& R2,
+				       Cell *& best)
+  {
+    CoordType d2 = 0;
+    int axis = level % N;
+    Node *other, *go;
+
+    if (x[axis] < node->value->coord[axis])
+      {
+	// The best is potentially in 0.
+	go = node->children[0];
+	other = node->children[1];
+      }
+    else
+      {
+	// If not it is in 1.
+	go = node->children[1];
+	other = node->children[0];
+	if (go == 0)
+	  {
+	    go = other;
+	    other = 0;
+	  }
+      }
+
+    if (go != 0)
+      {
+	recursiveNearest(go, level+1,
+			 x, R2,best);
+      }
+    else
+      {
+	CoordType thisR2 = computeDistance(node->value, x);
+	if (thisR2 < R2)
+	  {
+	    R2 = thisR2;
+	    best = node->value;
+	  }
+	return;
+      }
+    
+    // Check if current node is not the nearest
+    CoordType thisR2 = 
+      computeDistance(node->value, x);
+    
+    if (thisR2 < R2)
+      {
+	R2 = thisR2;
+	best = node->value;
+      }    
+    
+    // Now we found the best. We check whether the hypersphere
+    // intersect the hyperplane of the other branch
+
+    CoordType delta1;
+
+    delta1 = x[axis]-node->value->coord[axis];
+    if (delta1*delta1 < R2)
+      {
+	// The hypersphere intersects the hyperplane. Try the
+	// other branch
+	if (other != 0)
+	  {
+	    recursiveNearest(other, level+1, x, R2, best);
+	  }
+      }
+  }
+
+ template<int N, typename ValType, typename CType, typename CellSplitter>
+ KDCell<N,ValType,CType> *
+ KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbour(const coords& x)
+ {
+   CoordType R2 = INFINITY;
+   Cell *best = 0;
+
+   recursiveNearest(root, 0, x, R2, best);   
+   
+   return best;
+ }
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  void
+  KDTree<N,ValType,CType,CellSplitter>::recursiveMultipleNearest(RecursionMultipleInfo<N,ValType,CType>& info, Node *node,
+					       int level)
+  {
+    CoordType d2 = 0;
+    int axis = level % N;
+    Node *other, *go;
+
+    if (info.x[axis] < node->value->coord[axis])
+      {
+	// The best is potentially in 0.
+	go = node->children[0];
+	other = node->children[1];
+      }
+    else
+      {
+	// If not it is in 1.
+	go = node->children[1];
+	other = node->children[0];
+	//	if (go == 0)
+	//	  {
+	//	go = other;
+	//other = 0;
+	//}
+      }
+
+    if (go != 0)
+      {
+	recursiveMultipleNearest(info, go, level+1);
+      }
+    
+    // Check if current node is not the nearest
+    CoordType thisR2 = 
+      computeDistance(node->value, info.x);
+    info.queue.push(node->value, thisR2);
+    info.traversed++;
+    //    if (go == 0)
+    //      return;
+    
+    // Now we found the best. We check whether the hypersphere
+    // intersect the hyperplane of the other branch
+
+    CoordType delta1;
+
+    delta1 = info.x[axis]-node->value->coord[axis];
+    if (delta1*delta1 < info.queue.getMaxPriority())
+      {
+	// The hypersphere intersects the hyperplane. Try the
+	// other branch
+	if (other != 0)
+	  {
+	    recursiveMultipleNearest(info, other, level+1);
+	  }
+      }    
+  }
+
+ template<int N, typename ValType, typename CType, typename CellSplitter>
+ void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint32_t N2,
+						    Cell **cells)
+ {
+   RecursionMultipleInfo<N,ValType> info(x, cells, N2);
+   
+   for (int i = 0; i < N2; i++)
+     cells[i] = 0;
+
+   recursiveMultipleNearest(info, root, 0);
+
+   //   std::cout << "Traversed = " << info.traversed << std::endl;
+ }
+
+ template<int N, typename ValType, typename CType, typename CellSplitter>
+ void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint32_t N2,
+						    Cell **cells,
+						    CoordType *distances)
+ {
+   RecursionMultipleInfo<N,ValType> info(x, cells, N2);
+   
+   for (int i = 0; i < N2; i++)
+     cells[i] = 0;
+
+   recursiveMultipleNearest(info, root, 0);
+   memcpy(distances, info.queue.getPriorities(), sizeof(CoordType)*N2);
+
+   //   std::cout << "Traversed = " << info.traversed << std::endl;
+ }
+
+#ifdef __KD_TREE_SAVE_ON_DISK
+#define KDTREE_DISK_SIGNATURE "KDTREE"
+#define KDTREE_DISK_SIGNATURE_LEN 7
+  
+  template<int N, typename CType>
+  struct KDTreeOnDisk
+  {
+    long cell_id;
+    long children_node[2];
+    typename KDDef<N, CType>::KDCoordinates minBound, maxBound;
+  };  
+
+  struct KDTreeHeader
+  {
+    char id[KDTREE_DISK_SIGNATURE_LEN];
+    long nodesUsed, numCells;
+    long rootId;
+  };
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  void KDTree<N,ValType,CType,CellSplitter>::saveTree(std::ostream& o) const
+  {
+    KDTreeHeader h;
+
+    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++)
+      {
+	KDTreeOnDisk<N,CType> node_on_disk;
+
+	node_on_disk.cell_id = nodes[i].value - base_cell;
+	if (nodes[i].children[0] == 0)
+	  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] = -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));
+
+	o.write((char *)&node_on_disk, sizeof(node_on_disk));
+      }
+  }
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  KDTree<N,ValType,CType,CellSplitter>::KDTree(std::istream& in, Cell *cells, uint32_t Ncells)
+    throw (InvalidOnDiskKDTree)
+  {
+    KDTreeHeader h;
+
+    if (!in)
+      throw InvalidOnDiskKDTree();
+
+    in.read((char *)&h, sizeof(h));
+    if (!in || strncmp(h.id, KDTREE_DISK_SIGNATURE, KDTREE_DISK_SIGNATURE_LEN) != 0)
+      {
+         std::cerr << "KDTree Signature invalid" << std::endl;
+         throw InvalidOnDiskKDTree();
+      } 
+
+    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];
+    lastNode = h.nodesUsed;
+    numNodes = Ncells;
+
+    for (long i = 0; i < lastNode; i++)
+      {
+	KDTreeOnDisk<N,CType> node_on_disk;
+
+	in.read((char *)&node_on_disk, sizeof(node_on_disk));
+
+	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();
+	  }
+
+	nodes[i].value = base_cell + node_on_disk.cell_id;
+	if (node_on_disk.children_node[0] == -1)
+	  nodes[i].children[0] = 0;
+	else
+	  nodes[i].children[0] = nodes + node_on_disk.children_node[0];
+
+	if (node_on_disk.children_node[1] == -1)
+	  nodes[i].children[1] = 0;
+	else
+	  nodes[i].children[1] = nodes + node_on_disk.children_node[1];
+
+	memcpy(nodes[i].minBound, node_on_disk.minBound, sizeof(coords));
+	memcpy(nodes[i].maxBound, node_on_disk.maxBound, sizeof(coords));
+
+	int c;
+	for (c = 0; c < N; c++)
+	  if (nodes[i].value->coord[c] < nodes[i].minBound[c] ||
+	      nodes[i].value->coord[c] > nodes[i].maxBound[c])
+	    break;
+	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];    
+  }
+#endif
+
+};

external/cosmotool/src/newton.hpp

@@ -0,0 +1,30 @@
+#ifndef _COSMOTOOL_NEWTON_HPP
+#define _COSMOTOOL_NEWTON_HPP
+
+#include <cmath>
+
+namespace CosmoTool
+{
+  template<typename T, typename FunT>
+  T newtonSolver(T x0, FunT function, double residual = 1e-3)
+  {
+    T x, xold = x0;
+    T f_x = function.eval(x0);
+    T df_x = function.derivative(x0);
+    
+    x = xold - f_x/df_x;
+
+    while (std::abs(xold-x) > residual)
+      {	
+	xold = x;
+	f_x = function.eval(x);
+	df_x = function.derivative(x);
+	x = xold - f_x/df_x;
+      }
+	   
+    return x;
+  }
+
+};
+
+#endif

external/cosmotool/src/octTree.cpp

@@ -0,0 +1,175 @@
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "config.hpp"
+#include "octTree.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+//#define VERBOSE
+
+static uint32_t mypow(uint32_t i, uint32_t p)
+{
+  if (p == 0)
+    return 1;
+  else if (p == 1)
+    return i;
+
+  uint32_t k = p/2;
+  uint32_t j = mypow(i, k);
+  if (2*k==p)
+    return j*j;
+  else
+    return j*j*i;  
+}
+
+OctTree::OctTree(const FCoordinates *particles, octPtr numParticles, 
+		 uint32_t maxMeanTreeDepth, uint32_t maxAbsoluteDepth,
+		 uint32_t threshold)
+{
+  cout << "MeanTree=" << maxMeanTreeDepth << endl;
+  numCells = mypow(8, maxMeanTreeDepth);
+  assert(numCells < invalidOctCell);
+  //#ifdef VERBOSE
+  cerr << "Allocating " << numCells << " octtree cells" << endl;
+  //#endif
+  
+  for (int j = 0; j < 3; j++)
+    xMin[j] = particles[0][j];
+
+  for (octPtr i = 1; i < numParticles; i++)
+    {
+      for (int j = 0; j < 3; j++)
+	{
+	  if (particles[i][j] < xMin[j])
+	    xMin[j] = particles[i][j];      
+	}
+    }
+  
+  lenNorm = 0;
+  for (octPtr i = 0; i < numParticles; i++)
+    {
+      for (int j = 0; j < 3; j++)
+	{
+	  float delta = particles[i][j]-xMin[j];
+	  if (delta > lenNorm)
+	    lenNorm = delta;
+	}
+    }
+  cout << xMin[0] << " " << xMin[1] << " " << xMin[2] << " lNorm=" << lenNorm << endl; 
+
+  cells = new OctCell[numCells];
+  Lbox = (float)(octCoordTypeNorm+1);
+
+  cells[0].numberLeaves = 0;
+  for (int i = 0; i < 8; i++)
+    cells[0].children[i] = emptyOctCell;
+
+  lastNode = 1;
+  this->particles = particles;
+  this->numParticles = numParticles;
+  buildTree(maxAbsoluteDepth);
+  //#ifdef VERBOSE
+  cerr << "Used " << lastNode << " cells" << endl;
+  //#endif
+}
+
+OctTree::~OctTree()
+{
+  delete cells;
+}
+
+void OctTree::buildTree(uint32_t maxAbsoluteDepth)
+{
+  for (octPtr i = 0; i < numParticles; i++)
+    {
+      OctCoords rootCenter = { octCoordCenter, octCoordCenter, octCoordCenter };
+      insertParticle(0, // root node
+		     rootCenter,
+		     octCoordCenter,
+		     i,
+		     maxAbsoluteDepth);
+    }
+}
+
+
+void OctTree::insertParticle(octPtr node, 
+			     const OctCoords& icoord,
+			     octCoordType halfNodeLength,
+			     octPtr particleId,
+			     uint32_t maxAbsoluteDepth)
+{
+  
+#ifdef VERBOSE
+  cout << "Entering " << node << " (" << icoord[0] << "," << icoord[1] << "," << icoord[2] << ")" << endl;
+#endif
+  int octPos = 0;
+  int ipos[3] = { 0,0,0};
+  octPtr newNode;
+  OctCoords newCoord;
+  
+  cells[node].numberLeaves++;
+  if (maxAbsoluteDepth == 0)
+    {
+      // All children must be invalid.
+      for (int i = 0 ; i < 8; i++)
+	cells[node].children[i] = invalidOctCell;
+      
+      return;
+    }
+  
+  for (int j = 0; j < 3; j++)
+    {
+      float treePos = (particles[particleId][j]-xMin[j])*Lbox/lenNorm;
+      if ((octPtr)(treePos) > icoord[j])
+	{
+	  octPos |= (1 << j);
+	  ipos[j] = 1;
+	}
+    }
+  
+  if (cells[node].children[octPos] == emptyOctCell)
+    {
+      // Put the particle there.
+      cells[node].children[octPos] = particleId | octParticleMarker;
+      return;
+    }
+  
+  // If it is a node, explores it.
+  if (!(cells[node].children[octPos] & octParticleMarker))
+    {
+      assert(halfNodeLength >= 2);
+      // Compute coordinates
+      for (int j = 0; j < 3; j++)
+	newCoord[j] = icoord[j]+(2*ipos[j]-1)*halfNodeLength/2;
+      insertParticle(cells[node].children[octPos], newCoord, halfNodeLength/2, 
+		     particleId, maxAbsoluteDepth-1);
+      return;
+    }
+  
+  // We have a particle there. 
+  // Make a new node and insert the old particle into this node.
+  // Insert the new particle into the node also
+  // Finally put the node in place
+  
+  newNode = lastNode++;
+  assert(lastNode != numCells);
+  
+  for (int j = 0; j < 8; j++)
+    cells[newNode].children[j] = emptyOctCell;
+  cells[newNode].numberLeaves = 0;
+  
+  // Compute coordinates
+  for (int j = 0; j < 3; j++)
+    newCoord[j] = icoord[j]+(2*ipos[j]-1)*halfNodeLength/2;
+  
+  octPtr oldPartId = cells[node].children[octPos] & octParticleMask;
+  
+  insertParticle(newNode, newCoord, halfNodeLength/2,
+		 oldPartId, maxAbsoluteDepth-1);
+  insertParticle(newNode, newCoord, halfNodeLength/2,
+		 particleId, maxAbsoluteDepth-1);  
+  cells[node].children[octPos] = newNode;
+}
+

external/cosmotool/src/octTree.hpp

@@ -0,0 +1,145 @@
+#ifndef __COSMOTOOL_AMR_HPP
+#define __COSMOTOOL_AMR_HPP
+
+#include "config.hpp"
+
+namespace CosmoTool
+{
+
+  typedef uint32_t octPtr;
+  typedef uint16_t octCoordType;
+  static const uint16_t octCoordTypeNorm = 0xffff;
+  static const uint16_t octCoordCenter = 0x8000;
+
+  // This is also the root cell, but this one
+  // is never referenced (really ??).
+  static const octPtr invalidOctCell = 0xffffffff;
+  static const octPtr emptyOctCell = 0;
+  static const octPtr octParticleMarker = 0x80000000;
+  static const octPtr octParticleMask = 0x7fffffff;
+
+  typedef octCoordType OctCoords[3];
+
+  struct OctCell
+  {
+    octPtr numberLeaves;
+    octPtr children[8];
+  };
+
+  class OctTree
+  {
+  public:
+    OctTree(const FCoordinates *particles, octPtr numParticles,
+	    uint32_t maxTreeDepth,  uint32_t maxAbsoluteDepth,
+	    uint32_t threshold = 1);
+    ~OctTree();
+
+    void buildTree(uint32_t maxAbsoluteDepth);
+    void insertParticle(octPtr node, 
+			const OctCoords& icoord,
+			octCoordType halfNodeLength,
+			octPtr particleId,
+			uint32_t maxAbsoluteDepth);
+
+
+    octPtr getNumberLeaves() const {
+      return cells[0].numberLeaves;
+    }
+
+    static bool unconditioned(const FCoordinates&, octPtr, float, bool)
+    {
+      return true;
+    }
+
+    template<typename FunT>
+    void walkTree(FunT f)
+    {
+      walkTree(f, unconditioned);
+    }
+
+    template<typename FunT, typename CondT>
+    void walkTree(FunT f, CondT condition)
+    {
+      OctCoords rootCenter = { octCoordCenter, octCoordCenter, octCoordCenter };
+
+      walkTreeElements(f, condition, 0, rootCenter, octCoordCenter);
+    }
+
+    
+
+  protected:
+    const FCoordinates *particles;
+    octPtr numParticles;
+    OctCell *cells;
+    float Lbox;
+    octPtr lastNode;
+    octPtr numCells;  
+    float lenNorm;
+    float xMin[3];
+
+
+    static bool unconditioned()
+    {
+      return true;
+    }
+
+    template<typename FunT, typename CondT>
+    void walkTreeElements(FunT f, CondT condition,
+			  octPtr node, 
+			  const OctCoords& icoord,
+			  octCoordType halfNodeLength)
+    {
+      OctCoords newCoord;
+      FCoordinates center, realCenter;
+
+      for (int j = 0; j < 3; j++)
+	{
+	  center[j] = icoord[j]/(2.*octCoordCenter);
+	  realCenter[j] = xMin[j] + center[j]*lenNorm;
+	}
+
+      f(realCenter, cells[node].numberLeaves, lenNorm*halfNodeLength/(float)octCoordCenter,
+	cells[node].children[0] == invalidOctCell, // True if this is a meta-node
+	false);
+     
+      if (!condition(realCenter, cells[node].numberLeaves,
+		     lenNorm*halfNodeLength/(float)octCoordCenter,
+		     cells[node].children[0] == invalidOctCell))
+	return;
+       
+      for (int i = 0; i < 8; i++)
+	{
+	  octPtr newNode = cells[node].children[i];
+	  int ipos[3] = { (i&1), (i&2)>>1, (i&4)>>2 };
+	  
+	  if (newNode == emptyOctCell || newNode == invalidOctCell)
+	    continue;
+	  	  
+	  for (int j = 0; j < 3; j++)
+	    newCoord[j] = icoord[j]+(2*ipos[j]-1)*halfNodeLength/2;
+	  
+	  if (newNode & octParticleMarker)
+	    {
+	      for (int j = 0; j < 3; j++)
+		{
+		  center[j] = newCoord[j]/(2.*octCoordCenter);
+		  realCenter[j] = xMin[j] + lenNorm*center[j];
+		}
+
+	      f(realCenter, 
+		1, lenNorm*halfNodeLength/(2.*octCoordCenter),
+		false, true);
+	      continue;
+	    }
+
+	  walkTreeElements(f, condition, cells[node].children[i], newCoord, halfNodeLength/2);
+	}
+      
+    }
+
+  };
+  
+};
+
+
+#endif

external/cosmotool/src/pool.hpp

@@ -0,0 +1,166 @@
+#ifndef __COSMO_POOL_HPP
+#define __COSMO_POOL_HPP
+
+#include <list>
+#include "config.hpp"
+
+namespace CosmoTool
+{
+
+  template<typename T>
+  struct PoolNode
+  {
+    T *data;
+    uint32_t last_free, size;
+    PoolNode<T> *next;
+  };
+
+  template<typename T> class MemoryPool;
+
+  template<typename T>
+  class MemoryIterator
+  {
+  private:
+    friend  class MemoryPool<T>;
+
+    PoolNode<T> *cur, *previous;
+    uint32_t in_node;
+
+    MemoryIterator(PoolNode<T> *h) 
+    {
+      cur = h; 
+      previous = h; 
+      in_node = 0;
+    }
+  public:
+    MemoryIterator() { cur = 0; }
+    ~MemoryIterator() {}
+
+    const MemoryIterator& operator=(const MemoryIterator& i)
+    {
+      cur = i.cur;
+      previous = i.previous;
+      in_node = i.in_node;
+    }
+    
+    bool operator==(const MemoryIterator& i) const
+    {
+      return (cur == i.cur) && (in_node == i.in_node);
+    }
+
+    MemoryIterator& operator++()
+    {
+      if (cur == 0)
+	return *this;
+
+      in_node++;
+      if (in_node == cur->size)
+	{
+	  in_node = 0;
+	  previous = cur;
+	  cur = cur->next;
+	}
+      return *this;
+    }
+
+    T& operator*()
+    {
+      return cur->data[in_node];
+    }
+
+    T& operator->()
+    {
+      return cur->data[in_node];
+    }
+
+    
+    
+  };
+
+  // This is bare simple memory pools
+  template<typename T>
+  class MemoryPool
+  {
+  private:
+    uint32_t m_allocSize;
+    PoolNode<T> *head, *current;
+    typedef MemoryIterator<T> iterator;
+  public:
+    MemoryPool(uint32_t allocSize)
+      : m_allocSize(allocSize), head(0), current(0) {}
+
+    ~MemoryPool()
+    {
+      free_all();
+    }
+
+    void free_all()
+    {
+      PoolNode<T> *node = head;
+
+      while (node != 0)
+	{
+	  PoolNode<T> *next = node->next;
+	  	  
+	  delete[] node->data;
+	  delete node;
+	  node = next;
+	}
+      current = head = 0;
+    }
+
+    T *alloc()
+    {
+      T *ret = alloc_in_node();
+      return (ret == 0) ? alloc_new_in_node() : ret;
+    }
+    
+    iterator begin()
+    {
+      return iterator(head);
+    }
+    
+    iterator end()
+    {
+      return iterator(0);
+    }    
+
+  protected:
+    T *alloc_in_node()
+    {
+      if (current == 0 || current->last_free == current->size)
+	return 0;
+      return &current->data[current->last_free++];
+    }
+
+    T *alloc_new_in_node()
+    {
+      PoolNode<T> *newNode = new PoolNode<T>;
+      if (newNode == 0)
+	return 0;
+
+      newNode->last_free = 1;
+      newNode->size = m_allocSize;
+      newNode->data = new T[m_allocSize];
+      if (newNode->data == 0)
+	{
+	  delete newNode;
+	  return 0;
+	}
+      newNode->next = 0;
+
+      if (current == 0)
+	current = head = newNode;
+      else
+	{
+	  current->next = newNode;
+	  current = newNode;
+	}
+      return &newNode->data[0];
+    }
+
+  };
+  
+};
+
+#endif

external/cosmotool/src/powerSpectrum.cpp

@@ -0,0 +1,654 @@
+#include <vector>
+#include <algorithm>
+#include <iostream>
+#include <cmath>
+#include <fstream>
+#include <gsl/gsl_integration.h>
+#include "powerSpectrum.hpp"
+
+using namespace std;
+
+#define USE_GSL
+#define TOLERANCE 1e-6
+#define NUM_ITERATION 8000
+
+#define POWER_EFSTATHIOU 1
+#define HU_WIGGLES 2
+#define HU_BARYON 3
+#define OLD_POWERSPECTRUM 4
+#define POWER_BARDEEN 5
+#define POWER_SUGIYAMA 6
+#define POWER_BDM 7
+#define POWER_TEST 8
+
+#define POWER_SPECTRUM HU_WIGGLES
+
+namespace Cosmology {
+
+double n = 1.0;
+double K0 = 1;
+double V0 = 627;
+
+  double CMB_VECTOR[3] = {
+    56.759,
+    -540.02,
+    313.50
+  };
+
+
+// WMAP5
+double h = 0.719;
+double SIGMA8 = 0.77;
+double OMEGA_B = 0.043969;
+double OMEGA_C = 0.21259;
+
+
+// WMAP5-modification
+//double h = 0.719;
+//double SIGMA8 = 0.77;
+//double OMEGA_B = 0;
+//double OMEGA_C = 0.21259+0.043969;
+
+// LCDM STRAUSS ?
+//double h = 0.67;
+//double SIGMA8 = 0.67;
+//double OMEGA_B = 0;
+//double OMEGA_C = 0.30; 
+
+// SCDM STRAUSS
+//double h = 0.5;
+//double SIGMA8= 1.05;
+//double OMEGA_B = 0;
+//double OMEGA_C = 1;
+
+// Sugiyama test
+//double h = 0.5;
+//double SIGMA8= 0.5;//1.05;
+//double OMEGA_B = 0.0125*4;
+//double OMEGA_C = 0.1-OMEGA_B;
+
+// HU TEST
+//double h = 0.5;
+//double SIGMA8 = 0.5;
+//double OMEGA_B = 0.09;
+//double OMEGA_C = 0.21;
+
+// HDM STRAUSS
+//double h = 0.5;
+//double SIGMA8 = 0.86;
+//double OMEGA_B = 0;
+//double OMEGA_C = 1;
+
+// FOR "BEST FIT"
+//double h = 0.82;
+//double SIGMA8 = 0.76;
+//double OMEGA_B = 0.043969;
+//double OMEGA_C = 0.15259;
+
+
+// FOR JUSZKIEWICZ CHECKING (CDM) ! WARNING ! He smoothes 
+// with a gaussian filter the density field, i.e. one has
+// to multiply P(k) by exp(-k^2 R^2) with R the radius
+// of the filter. Dammit !
+//double h = 0.5;
+//double SIGMA8=1/2.5;
+//double OMEGA_B=0.;
+//double OMEGA_C=1;
+//#define JUSZKIEWICZ_PATCH
+//#define RJUSZ 6.0
+
+// (BDM)
+//double h = 0.5;
+//double SIGMA8=1;
+//double OMEGA_B=0.0;
+//double OMEGA_C=0.4;
+
+
+// FOR HU CHECKING
+//double h = 0.5;
+//double SIGMA8= 1;
+//double OMEGA_B=0.09;
+//double OMEGA_C=0.21;
+
+
+double OMEGA_0 = OMEGA_B+OMEGA_C;
+double Omega = OMEGA_0;
+double Theta_27 = 2.728 / 2.7;
+double beta = pow(OMEGA_0, 5./9);
+double OmegaEff = OMEGA_0 * h * h;
+double Gamma0 = OMEGA_0 * h * h;
+
+
+/*
+ * This is \hat{tophat}
+ */
+double tophatFilter(double u)
+{
+  if (u != 0)    
+    return 3 / (u*u*u) * (sin(u) - u * cos(u));
+  else
+    return 1;
+}
+
+/*
+ * This is \tilde{w}
+ */
+double externalFilter(double u)
+{
+  if (u != 0)
+    return 1 - sin(u)/u;
+  return 0.;
+}
+
+double powC(double q, double alpha_c)
+{
+  return 14.2 / alpha_c + 386 / (1 + 69.9 * pow(q, 1.08));
+}
+
+double T_tilde_0(double q, double alpha_c, double beta_c)
+{
+  double a = log(M_E + 1.8 * beta_c * q);
+  return a / ( a  + powC(q, alpha_c) * q * q);
+}
+
+double j_0(double x)
+{
+  if (x == 0)
+    return 1.0;
+  return sin(x)/x;
+}
+
+double powG(double y)
+{
+  double a = sqrt(1 + y);
+
+  return y * (-6 * a + (2 + 3 * y) *log((a + 1)/(a - 1)));
+}
+
+/*
+ * This function returns the power spectrum evaluated at k (in Mpc, not in Mpc/h).
+ */
+double powerSpectrum(double k, double normPower)
+{
+#if POWER_SPECTRUM == POWER_EFSTATHIOU
+  double a = 6.4/Gamma0;
+  double b = 3/Gamma0;
+  double c = 1.7/Gamma0;
+  double nu = 1.13;
+  
+  double f = (a*k) + pow(b*k,1.5) + pow(c*k,2);
+
+  // EFSTATHIOU  ET AL. (1992)
+  return normPower * pow(k,n) * pow(1+pow(f,nu),(-2/nu));
+#endif
+
+  // EISENSTEIN ET HU (1998)
+  // FULL POWER SPECTRUM WITH BARYONS AND WIGGLES
+
+#if POWER_SPECTRUM == HU_WIGGLES
+  // EISENSTEIN ET HU (1998)
+  // FULL POWER SPECTRUM WITH BARYONS AND WIGGLES
+
+  double k_silk = 1.6 * pow(OMEGA_B * h * h, 0.52) * pow(OmegaEff, 0.73) * (1 + pow(10.4 * OmegaEff, -0.95));
+  double z_eq = 2.50e4 * OmegaEff * pow(Theta_27, -4);
+  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double f = 1 / (1 + pow(k * s / 5.4, 4));
+  double k_eq = 7.46e-2 * OmegaEff * pow(Theta_27, -2);
+  double a1 = pow(46.9 * OmegaEff, 0.670) * (1 + pow(32.1 * OmegaEff, -0.532));
+  double a2 = pow(12.0 * OmegaEff, 0.424) * (1 + pow(45.0 * OmegaEff, -0.582));
+  double alpha_c = pow(a1, -OMEGA_B/ OMEGA_0) * pow(a2, -pow(OMEGA_B / OMEGA_0, 3));  
+
+  double q = k / (13.41 * k_eq);
+  double b1_betac = 0.944 * 1/(1 + pow(458 * OmegaEff, -0.708));
+  double b2_betac = pow(0.395 * OmegaEff, -0.0266);
+  double beta_c = 1/ ( 1 + b1_betac * (pow(OMEGA_C / OMEGA_0, b2_betac) - 1)   );
+  double T_c = f * T_tilde_0(q, 1, beta_c) + (1 - f) * T_tilde_0(q, alpha_c, beta_c);
+
+  double b1_zd = 0.313 * pow(OmegaEff, -0.419) * (1 + 0.607 * pow(OmegaEff, 0.674));
+  double b2_zd = 0.238 * pow(OmegaEff, 0.223);
+  double z_d = 1291 * pow(OmegaEff, 0.251) / (1 + 0.659 * pow(OmegaEff, 0.828)) * (1 + b1_zd * pow(OmegaEff, b2_zd));
+  double R_d = 31.5 * OMEGA_B * h * h * pow(Theta_27, -4) * 1e3 / z_d;
+
+  double alpha_b = 2.07 * k_eq * s * pow(1 + R_d, -0.75) * powG((1 + z_eq)/(1 + z_d));
+  double beta_b = 0.5 + OMEGA_B / OMEGA_0 + (3 - 2 * OMEGA_B / OMEGA_0) * sqrt(pow(17.2 * OmegaEff, 2) + 1);
+  double beta_node = 8.41 * pow(OmegaEff, 0.435);
+  double s_tilde = s * pow(1 + pow(beta_node / (k * s), 3), -1./3);
+
+  double T_b = (T_tilde_0(q, 1, 1) / (1 + pow(k * s / 5.2, 2)) + alpha_b / (1 + pow(beta_b / (k * s), 3)) * exp(-pow(k/k_silk, 1.4))) * j_0(k * s_tilde);  
+
+  double T_k = OMEGA_B/OMEGA_0 * T_b + OMEGA_C/OMEGA_0 * T_c;  
+
+  return normPower * pow(k,n) * T_k * T_k;
+#endif
+
+  // EISENSTEIN ET AL. (2008), SHAPED POWER SPECTRUM WITH BARYON, WITHOUT WIGGLES
+#if POWER_SPECTRUM == HU_BARYON
+  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2);
+  double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4)));
+  double q = k/(h*GammaEff) * pow(Theta_27, 2);
+  double L_0 = log(2 * M_E + 1.8 * q);
+  double C_0 = 14.2 + 731 / (1 + 62.5 * q);
+  double T0 = L_0 / (L_0 + C_0 * q * q);
+  
+  return normPower * pow(k,n) * T0 * T0;
+#endif
+
+#if POWER_SPECTRUM == OLD_POWERSPECTRUM
+  // OLD FUNCTION: 
+  static const double l = 1/(Omega * h*h);
+  static const double alpha = 1.7 * l, beta = 9.0 * pow(l, 1.5), gamma = l*l;
+  return normPower * pow(k,n) * pow(1 + alpha * k +  beta * pow(k,1.5) + gamma *k*k,-2);
+#endif
+
+#if POWER_SPECTRUM == POWER_SUGIYAMA
+  double q = k * Theta_27*Theta_27 / (OmegaEff * exp(-OMEGA_B - sqrt(h/0.5)*OMEGA_B/OMEGA_0));
+  double L0 = log(2*M_E + 1.8 * q);
+  double C0 = 14.2 + 731 / (1 + 62.5 * q);
+  double T_k = L0 / (L0 + C0 * q*q);
+
+      //  double poly = 1 + 3.89 * q + pow(16.1*q,2) + pow(5.46*q,3) + pow(6.71*q,4);
+      //  double T_k = log(1+2.34*q)/(2.34*q) * pow(poly,-0.25);
+
+  return normPower * pow(k,n) * T_k * T_k;
+
+#endif
+
+#if POWER_SPECTRUM == POWER_BARDEEN
+  double q = k / (OmegaEff);
+
+  double poly = 1 + 3.89 * q + pow(16.1*q,2) + pow(5.46*q,3) + pow(6.71*q,4);
+  double T_k = log(1+2.34*q)/(2.34*q) * pow(poly,-0.25);
+
+  return normPower * pow(k,n) * T_k * T_k;
+
+#endif
+
+#if POWER_SPECTRUM == POWER_BDM
+
+  k /= h*h;
+  double alpha1 = 190;
+  double Gmu = 4.6;
+  double alpha2 = 11.5;
+  double alpha3 = 11;
+  double alpha4 = 12.55;
+  double alpha5 = 0.0004;
+  return normPower*k*alpha1*alpha1*Gmu*Gmu/(1+(alpha2*k)+pow(alpha3*k,2)+pow(alpha4*k,3))*pow(1+pow(alpha5/k,2), -2);
+#endif
+
+#if POWER_SPECTRUM == POWER_TEST
+  return 1/(1+k*k);
+#endif
+}
+
+/*
+ * This function computes the normalization of the power spectrum. It requests
+ * a sigma8 (density fluctuations within 8 Mpc/h)
+ */
+double gslPowSpecNorm(double k, void *params)
+{
+  double f = tophatFilter(k*8.0/h);
+  return powerSpectrum(k, 1.0)*k*k*f*f;
+}
+
+double computePowSpecNorm(double sigma8)
+{
+  int Nsteps = 30000;
+  double normVal = 0;
+
+#ifndef USE_GSL
+  for (int i = 1; i <= Nsteps; i++)
+    {
+      double t = i * 1.0/(Nsteps+1);
+      // Change of variable !
+      double k = (1-t)/t * K0;
+
+      // The filter
+      double filter_val = tophatFilter(k*8.0/h);
+      // The powerspectrum
+      double powVal = powerSpectrum(k, 1.0);
+      // Multiply by the tophat filter
+      powVal *= filter_val*filter_val;
+     
+      powVal *= k*k;
+      
+      // Account for change of variable
+      powVal /= (t*t);
+  
+      // Integrate !
+      normVal += powVal;
+    }
+
+  normVal /= 2*M_PI*M_PI;
+
+  // The dt element
+  normVal *= 1.0/(Nsteps+1) * K0;
+#else
+  double abserr;
+  gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION);
+  gsl_function f;
+
+  f.function = gslPowSpecNorm;
+  f.params = 0;
+
+  gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr);
+
+  gsl_integration_workspace_free(w);
+
+  normVal /= (2*M_PI*M_PI);
+#endif
+
+  return sigma8*sigma8/normVal;
+}
+
+/*
+ * This function computes the variance of the Local Group velocity components
+ * for a survey which depth is topHatRad1 (in Mpc/h). This variance should
+ * be multiplied by (H \beta)^2 to be equal to a velocity^2.
+ */
+double gslVariance(double k, void *params)
+{
+  double R1 = *(double *)params;
+  double f = externalFilter(k * R1 / h);
+
+  double a = f*f;
+
+#ifdef JUSZKIEWICZ_PATCH
+  a *= exp(-k*k*(RJUSZ*RJUSZ/(h*h)));
+#endif
+  a *= powerSpectrum(k, 1.0);
+
+  return a;
+}
+
+double computeVariance(double powNorm,  double topHatRad1)
+{
+  int Nsteps = 100000;
+  double varVal = 0;
+
+#ifndef USE_GSL
+  
+  for (int i = 1; i <= Nsteps; i++)
+    {
+      double t = i * 1.0/(Nsteps+1);
+      // Change of variable !
+      double k = (1-t)/t * K0;
+      double powVal = powerSpectrum(k, powNorm);
+
+      double filter1Val = externalFilter(k*topHatRad1/h);
+
+#ifdef JUSZKIEWICZ_PATCH
+      powVal *= exp(-k*k*(RJUSZ*RJUSZ/(h*h)));
+#endif
+      powVal *= filter1Val*filter1Val;
+      
+      powVal /= (t*t);
+
+      varVal += powVal;
+    }
+
+  varVal *= 1.0/(Nsteps) * K0;
+  varVal *= 1.0/(6*M_PI*M_PI);
+
+#else
+
+  double abserr;
+  gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION);
+  gsl_function f;
+  
+  f.function = gslVariance;
+  f.params = &topHatRad1;
+
+  gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &varVal, &abserr);
+
+  gsl_integration_workspace_free(w);
+
+  varVal *= powNorm/(6*M_PI*M_PI);
+
+#endif
+
+  return varVal;
+}
+
+/*
+ * This function computes the same quantity as computeVariance but
+ * for a survey infinitely deep.
+ */
+double gslVarianceZero(double k, void *params)
+{
+  double a = 1.0;
+
+#ifdef JUSZKIEWICZ_PATCH
+  a *= exp(-k*k*(RJUSZ*RJUSZ/(h*h)));
+#endif
+  a *= powerSpectrum(k, 1.0);
+
+  return a;
+}
+
+double computeVarianceZero(double powNorm)
+{
+  int Nsteps = 100000;
+  double varVal = 0;
+  
+#ifndef USE_GSL
+
+  for (int i = 1; i <= Nsteps; i++)
+    {
+      double t = i * 1.0/(Nsteps+1);
+      // Change of variable !
+      double k = (1-t)/t * K0;
+      double powVal = powerSpectrum(k, powNorm);
+      
+#ifdef JUSZKIEWICZ_PATCH
+      powVal *= exp(-k*k*(RJUSZ*RJUSZ/h*h));
+#endif
+
+      powVal /= (t*t);
+
+      varVal += powVal;
+    }
+
+  varVal *= 1.0/(Nsteps+1) * K0;
+  varVal *= 1.0/(6*M_PI*M_PI);
+
+#else
+ 
+  double abserr;
+  gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION);
+  gsl_function f;
+  
+  f.function = gslVarianceZero;
+  f.params = 0;
+
+  gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &varVal, &abserr);
+
+  gsl_integration_workspace_free(w);
+
+  varVal *= powNorm/(6*M_PI*M_PI);
+#endif
+
+  return varVal;
+}
+
+/*
+ * This function computes the correlation between the infinitely deep
+ * velocity of the Local Group and the one estimated from a survey
+ * which depth is topHatRad1.
+ * This corresponds to \gamma.
+ * This quantity must be multiplied by H \beta to be equal to a velocity^2.
+ */
+double gslCorrel(double k, void *params)
+{
+  double R1 = ((double *)params)[0];
+  double a = externalFilter(k * R1 / h);// * externalFilter(k * R2 / h);
+
+#ifdef JUSZKIEWICZ_PATCH
+  a *= exp(-k*k*(RJUSZ*RJUSZ/(h*h)));
+#endif
+  a *= powerSpectrum(k, 1.0);
+
+  return a;
+}
+
+double gslCorrelBis(double t, void *params)
+{
+  double k = (1-t)/t;
+  double v = gslCorrel(k, params);
+
+  return v/(t*t);
+}
+
+double gslCorrel2(double k, void *params)
+{
+  double R1 = ((double *)params)[0];
+  double R2 = ((double *)params)[1];
+  double a = externalFilter(k * R1 / h) * externalFilter(k * R2 / h);        
+  
+#ifdef JUSZKIEWICZ_PATCH
+  a *= exp(-k*k*(RJUSZ*RJUSZ/(h*h)));
+#endif
+  a *= powerSpectrum(k, 1.0) ;
+  
+  return a;
+}
+
+double gslCorrel2bis(double t, void *params)
+{
+  double k = (1-t)/t;
+  double v = gslCorrel2(k, params);
+
+  return v/(t*t);
+}
+
+
+double computeCorrel(double powNorm,  double topHatRad1)
+{
+  int Nsteps = 100000;
+  double varVal = 0;
+  
+#ifndef USE_GSL
+
+  for (int i = 1; i <= Nsteps; i++)
+    {
+      double t = i * 1.0/(Nsteps+1);
+      // Change of variable !
+      double k = (1-t)/t * K0;
+      double powVal = powerSpectrum(k, powNorm);
+
+      double filter1Val = externalFilter(k*topHatRad1/h);
+
+#ifdef JUSZKIEWICZ_PATCH
+      powVal*=exp(-k*k*(RJUSZ*RJUSZ/h*h));
+#endif
+
+      powVal *= filter1Val;
+      
+      powVal /= (t*t);
+
+      varVal += powVal;
+    }
+
+  varVal *= 1.0/(Nsteps) * K0;
+  varVal *= 1.0/(6*M_PI*M_PI);
+#else
+
+  double abserr;
+  gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION);
+  gsl_function f;
+  
+  f.params = &topHatRad1;
+
+#if 1
+  f.function = gslCorrelBis;
+  gsl_integration_qag (&f, 0, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &varVal, &abserr);
+#else
+  f.function = gslCorrel;
+  gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &varVal, &abserr);
+#endif
+
+  gsl_integration_workspace_free(w);
+
+  varVal *= powNorm/(6*M_PI*M_PI);
+
+#endif
+  return varVal;
+}
+
+/*
+ * This function computes the correlation between the infinitely deep
+ * velocity of the Local Group and the one estimated from a survey
+ * which depth is topHatRad1.
+ * This corresponds to \gamma.
+ * This quantity must be multiplied by H \beta to be equal to a velocity^2.
+ */
+double computeCorrel2(double powNorm,  double topHatRad1, double topHatRad2)
+{
+  int Nsteps = 100000;
+  double varVal = 0;
+#ifndef USE_GSL
+  
+  for (int i = 1; i <= Nsteps; i++)
+    {
+      double t = i * 1.0/(Nsteps+1);
+      // Change of variable !
+      double k = (1-t)/t * K0;
+      double powVal = powerSpectrum(k, powNorm);
+
+      double filter1Val = externalFilter(k*topHatRad1/h);
+      double filter2Val = externalFilter(k*topHatRad2/h);
+
+      powVal *= filter1Val * filter2Val;
+      
+      powVal /= (t*t);
+
+      varVal += powVal;
+    }
+
+  varVal *= 1.0/(Nsteps) * K0;
+  varVal *= 1.0/(6*M_PI*M_PI);
+
+#else
+
+  double abserr;
+  gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION);
+  gsl_function f;
+  double rads[] = {topHatRad1, topHatRad2 };
+
+  f.params = rads;
+
+#if 1  
+  f.function = gslCorrel2bis;
+  gsl_integration_qag (&f, 0, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &varVal, &abserr);
+#else
+  f.function = gslCorrel2;
+  gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &varVal, &abserr);
+#endif
+
+  gsl_integration_workspace_free(w);
+
+  varVal *= powNorm/(6*M_PI*M_PI);
+
+#endif
+
+  return varVal;
+}
+
+
+  void updateCosmology()
+  {
+    OMEGA_0 = OMEGA_B+OMEGA_C;
+    Omega = OMEGA_0;
+    Theta_27 = 2.728 / 2.7;
+    beta = pow(OMEGA_0, 5./9);
+    OmegaEff = OMEGA_0 * h * h;
+    Gamma0 = OMEGA_0 * h * h;
+
+#if 0
+    cout << "Cosmology is :" << endl
+	 << "  O0=" << OMEGA_0 << " Theta=" << Theta_27 << " beta=" << beta << " h=" << h << " G0=" << Gamma0 << endl
+	 << "  OmegaB=" << OMEGA_B << " Omega_C=" << OMEGA_C << endl;
+#endif
+  }
+
+};

external/cosmotool/src/powerSpectrum.hpp

@@ -0,0 +1,35 @@
+#ifndef _POWERSPECTRUM_HPP
+#define _POWERSPECTRUM_HPP
+
+namespace Cosmology {
+
+  extern double n;
+  extern double K0;
+  extern double V0;
+
+  extern double CMB_VECTOR[3];
+
+  // WMAP5
+  extern double h;
+  extern double SIGMA8;
+  extern double OMEGA_B;
+  extern double OMEGA_C;
+
+  extern double OMEGA_0;
+  extern double Omega;
+  extern double Theta_27;
+  extern double beta;
+  extern double OmegaEff;
+  extern double Gamma0;
+
+  void updateCosmology();
+  double tophatFilter(double u);
+  double powerSpectrum(double k, double normPower);
+  double computePowSpecNorm(double sigma8);
+  double computeVariance(double powNorm,  double topHatRad1);
+  double computeVarianceZero(double powNorm);
+  double computeCorrel(double powNorm,  double topHatRad1);
+  double computeCorrel2(double powNorm,  double topHatRad1, double topHatRad2);
+};
+
+#endif

external/cosmotool/src/sparseGrid.hpp

@@ -0,0 +1,87 @@
+#ifndef __VOID_SPARSE_GRID_HPP
+#define __VOID_SPARSE_GRID_HPP
+
+#include <cmath>
+#include "config.hpp"
+#include <list>
+
+namespace CosmoTool 
+{
+  
+  template <typename T>
+  struct bracketAccessor
+  {
+    typedef typename T::value_type result_type;
+
+    result_type
+    operator()(T const& V, size_t const N) const throw ()
+    {
+      return V[N];
+    }
+  };
+
+  template<typename T>
+  struct GridElement
+  {
+    std::list<T> eList;
+  };
+
+  template<typename T, typename CType = double>
+  struct GridDesc
+  {
+    GridElement<T> *grid;
+    uint32_t gridRes;
+    CType boxSize;
+  };
+
+  template<typename T, int N, typename CType = double>
+  class SparseIterator
+  {
+  protected:
+    typedef CType coord[N];
+
+    int32_t index_min[N], index_max[N], index_cur[N];
+    GridElement<T> *eCur;
+    GridDesc<T> desc;
+    typename std::list<T>::iterator eIterator;
+
+    GridElement<T> *getElement();
+  public:
+    SparseIterator(GridDesc<T>& d, coord& c, CType R);
+    SparseIterator();
+    ~SparseIterator();
+    
+    SparseIterator& operator++();
+
+    bool last();
+
+    T& operator*() { return *eIterator; }
+    T *operator->() { return eIterator.operator->(); }
+  };
+
+  template<typename T, int N, typename CType = double,
+	   typename CAccessor = bracketAccessor<T> >
+  class SparseGrid
+  {
+  public:
+    typedef CType coord[N];
+    typedef SparseIterator<T, N, CType> iterator;
+    CAccessor acc;
+
+    SparseGrid(CType L, uint32_t gridRes);
+    ~SparseGrid();
+    
+    void addElementInGrid(T& e);
+
+    iterator locateElements(coord& c, CType R);
+    void clear();
+
+  protected:
+    GridDesc<T> grid;
+  };
+
+};
+
+#include "sparseGrid.tcc"
+
+#endif

external/cosmotool/src/sparseGrid.tcc

@@ -0,0 +1,145 @@
+#include <cassert>
+#include <cstring>
+
+namespace CosmoTool {
+   
+   template<typename T, int N, typename CType,
+	    typename CAccessor>
+   SparseGrid<T,N,CType,CAccessor>::SparseGrid(CType L, uint32_t gridRes)
+   {
+     uint32_t Ntot = 1;
+     for (int i = 0; i < N; i++)
+       Ntot *= gridRes;
+
+     grid.grid = new GridElement<T>[Ntot];
+     grid.boxSize = L;
+     grid.gridRes = gridRes;
+   }
+
+   template<typename T, int N, typename CType,
+	    typename CAccessor>
+   SparseGrid<T,N,CType,CAccessor>::~SparseGrid()
+   {
+     delete[] grid.grid;
+   }
+    
+   template<typename T, int N, typename CType,
+	    typename CAccessor>
+   void SparseGrid<T,N,CType,CAccessor>::addElementInGrid(T& e)
+   {
+     uint32_t index[N];
+     uint32_t idx = 0;
+
+     for (int i = N-1; i >= 0; i--)
+       {
+	 CType x = acc(e, i) * grid.gridRes / grid.boxSize;
+	 
+	 index[i] = (int32_t)std::floor((double)x);
+	 assert(index[i] >= 0);
+	 assert(index[i] < grid.gridRes);
+	 idx = (idx*grid.gridRes) + index[i];
+       }
+     assert(idx < grid.gridRes*grid.gridRes*grid.gridRes);
+
+     grid.grid[idx].eList.push_back(e);
+   }
+
+   template<typename T, int N, typename CType,
+	    typename CAccessor>
+   SparseIterator<T,N,CType>
+      SparseGrid<T,N,CType,CAccessor>::locateElements(coord& c, CType R)
+   {
+     return SparseIterator<T,N,CType>(grid, c, R);
+   }
+
+  template<typename T, int N, typename CType,
+	   typename CAccessor>
+  void SparseGrid<T,N,CType,CAccessor>::clear()
+  {
+    uint32_t N3 = grid.gridRes*grid.gridRes*grid.gridRes;
+
+    for (uint32_t i = 0; i < N3; i++)
+      grid.grid[i].eList.clear();
+  }
+
+
+  template<typename T, int N, typename CType>
+  SparseIterator<T,N,CType>::SparseIterator(GridDesc<T>& d, coord& c, CType R)
+  {
+    desc = d;
+    
+    for (uint32_t i = 0; i < N; i++)
+      {
+	CType x_min = (c[i] - R) * d.gridRes / d.boxSize;
+	CType x_max = (c[i] + R) * d.gridRes / d.boxSize;
+	
+	index_min[i] = (int32_t) std::floor((double)x_min);
+	index_max[i] = (int32_t) std::ceil((double)x_max) + 1;
+	index_cur[i] = index_min[i];
+      }
+
+    eCur = getElement();
+    eIterator = eCur->eList.begin();
+    if (eIterator == eCur->eList.end())
+      operator++();
+  }
+  
+  template<typename T, int N, typename CType>
+  SparseIterator<T,N,CType>::~SparseIterator()
+  {
+  }
+  
+  template<typename T, int N, typename CType>
+  bool SparseIterator<T,N,CType>::last()
+  {
+    return (index_cur[N-1] == index_max[N-1]);
+  }
+
+ 
+  template<typename T, int N, typename CType>
+  GridElement<T> *SparseIterator<T,N,CType>::getElement()
+  {
+    uint32_t idx = 0;
+    for (int i = (N-1); i >= 0; i--)
+      {
+	int32_t k = (index_cur[i] + desc.gridRes) % desc.gridRes;	
+	idx = (idx*desc.gridRes) + k;
+      }
+    assert(idx < desc.gridRes*desc.gridRes*desc.gridRes);
+    
+    return &desc.grid[idx];
+  }
+
+  template<typename T, int N, typename CType>
+  SparseIterator<T,N,CType>& SparseIterator<T,N,CType>::operator++()
+  {
+    if (last())
+      return *this;
+    ++eIterator;
+    if (eIterator != eCur->eList.end())
+      return *this;
+      
+    do
+      {		
+	index_cur[0]++;
+	for (int i = 0; i < (N-1); i++)
+	  if (index_cur[i] == index_max[i])
+	    {
+	      index_cur[i] = index_min[i];
+	      index_cur[i+1]++;
+	    }
+	
+	if (last())
+	  return *this;
+	
+	eCur = getElement();
+	eIterator = eCur->eList.begin();
+      }
+    while (eIterator == eCur->eList.end());
+
+    return *this;
+  }
+ 
+
+
+};

external/cosmotool/src/sphSmooth.hpp

@@ -0,0 +1,105 @@
+#ifndef __COSMOTOOL_SPH_SMOOTH_HPP
+#define __COSMOTOOL_SPH_SMOOTH_HPP
+
+#include "config.hpp"
+#include "mykdtree.hpp"
+
+namespace CosmoTool
+{
+  template <typename ValType, int Ndims = NUMDIMS>
+  class SPHSmooth
+  {
+  public:
+    typedef struct
+    {
+      ComputePrecision weight;
+      ValType pValue;
+    } FullType;
+
+    typedef KDTree<Ndims,FullType> SPHTree;
+    typedef KDTreeNode<Ndims,FullType> SPHNode;
+    typedef KDCell<Ndims,FullType>  SPHCell;
+    typedef typename KDTree<Ndims,FullType>::CoordType CoordType;
+
+    typedef ComputePrecision (*computeParticleValue)(const ValType& t);
+    typedef void (*runParticleValue)(ValType& t);
+
+  public:
+    SPHSmooth(SPHTree *tree, uint32_t Nsph);
+    virtual ~SPHSmooth();   
+
+    void fetchNeighbours(const typename SPHTree::coords& c);
+    void fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNsph);
+    void fetchNeighboursOnVolume(const typename SPHTree::coords& c, ComputePrecision radius);
+    const typename SPHTree::coords& getCurrentCenter() const
+    {
+      return currentCenter;
+    }
+
+    ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c,
+					  computeParticleValue fun);    
+    ComputePrecision computeInterpolatedValue(const typename SPHTree::coords& c,
+					      computeParticleValue fun);    
+    ComputePrecision getMaxDistance(const typename SPHTree::coords& c, 
+				    SPHNode *node) const;
+
+    ComputePrecision getSmoothingLen() const
+    {
+      return smoothRadius;
+    }
+
+    // TO USE WITH EXTREME CARE !
+    void setSmoothingLen(ComputePrecision len)
+    {
+      smoothRadius = len;
+    }
+    // END
+
+    void runForEachNeighbour(runParticleValue fun);
+    void addGridSite(const typename SPHTree::coords& c);
+
+    bool hasNeighbours() const;
+
+    virtual ComputePrecision getKernel(ComputePrecision d) const;    
+
+    SPHTree *getTree() { return tree; }
+
+    void changeNgb(uint32_t newMax) {
+      delete[] ngb;
+      delete[] distances;
+      ngb = new SPHCell *[newMax];
+      distances = new CoordType[newMax];
+      maxNgb = newMax;
+    }
+
+    uint32_t getCurrent() const { return currentNgb; }
+
+    uint32_t getNgb() const { return maxNgb; }
+ 
+  protected:
+    SPHCell **ngb;
+    CoordType *distances;
+    uint32_t Nsph;
+    uint32_t deltaNsph;
+    uint32_t maxNgb;
+    uint32_t currentNgb;
+    SPHTree *tree;
+    ComputePrecision smoothRadius;
+    typename SPHTree::coords currentCenter;
+    
+    ComputePrecision computeWValue(const typename SPHTree::coords & c,
+				   SPHCell& cell,
+				   CoordType d,
+				   computeParticleValue fun);
+    void runUnrollNode(SPHNode *node,
+		       runParticleValue fun);
+  };
+
+  template<class ValType1, class ValType2, int Ndims>
+  bool operator<(const SPHSmooth<ValType1, Ndims>& s1, const SPHSmooth<ValType2, Ndims>& s2);
+
+};
+
+#include "sphSmooth.tcc"
+
+#endif

external/cosmotool/src/sphSmooth.tcc

@@ -0,0 +1,214 @@
+#include <cmath>
+
+namespace CosmoTool {
+
+template<typename ValType, int Ndims>
+SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph)
+{
+  this->Nsph = Nsph;
+  this->tree = tree;
+  this->currentNgb = 0;
+
+  this->maxNgb = Nsph;
+  ngb = new SPHCell *[maxNgb];  
+  distances = new CoordType[maxNgb];
+}
+
+template<typename ValType, int Ndims>
+SPHSmooth<ValType,Ndims>::~SPHSmooth()
+{
+  delete[] ngb;
+  delete[] distances;
+}
+
+template<typename ValType, int Ndims>
+ComputePrecision SPHSmooth<ValType,Ndims>::computeWValue(const typename SPHTree::coords& c,
+							 SPHCell& cell,
+							 CoordType d,
+							 computeParticleValue fun)
+{
+  CoordType weight;
+
+  d /= smoothRadius;
+  weight = getKernel(d);
+
+  if (cell.val.weight != 0)
+    return weight * fun(cell.val.pValue) / cell.val.weight;
+  else
+    return 0;
+}
+
+template<typename ValType, int Ndims>
+void
+SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNngb)
+{
+  ComputePrecision d2, max_dist = 0;
+  uint32_t requested = newNngb;
+
+  if (requested > maxNgb)
+    {
+      delete[] ngb;
+      delete[] distances;
+
+      maxNgb = requested;
+      ngb = new SPHCell *[maxNgb];
+      distances = new CoordType[maxNgb];
+    }
+
+  memcpy(currentCenter, c, sizeof(c));
+  tree->getNearestNeighbours(c, requested, ngb, distances);
+
+  currentNgb = 0;
+  for (uint32_t i = 0; i < requested && ngb[i] != 0; i++,currentNgb++)
+    {
+      distances[i] = sqrt(distances[i]);
+      d2 = distances[i];
+      if (d2 > max_dist)
+        max_dist = d2;
+    }
+
+  smoothRadius = max_dist / 2;  
+}
+
+template<typename ValType, int Ndims>
+void
+SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c)
+{
+  ComputePrecision d2, max_dist = 0;
+  uint32_t requested = Nsph;
+
+  memcpy(currentCenter, c, sizeof(c));
+  tree->getNearestNeighbours(c, requested, ngb, distances);
+
+  currentNgb = 0;
+  for (uint32_t i = 0; i < requested && ngb[i] != 0; i++,currentNgb++)
+    {
+      distances[i] = sqrt(distances[i]);
+      d2 = distances[i];
+      if (d2 > max_dist)
+	max_dist = d2;
+    }
+
+  smoothRadius = max_dist / 2;  
+}					
+
+template<typename ValType, int Ndims>
+void
+SPHSmooth<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords& c,
+						  ComputePrecision radius)
+{
+  uint32_t numPart;
+  ComputePrecision d2, max_dist = 0;
+
+  memcpy(currentCenter, c, sizeof(c));
+
+  currentNgb = tree->getIntersection(c, radius, ngb, distances,
+				  maxNgb);
+
+  for (uint32_t i = 0; i < currentNgb; i++)
+    {
+      distances[i] = sqrt(distances[i]);
+      d2 = distances[i];
+      if (d2 > max_dist)
+	max_dist = d2;
+    }
+  smoothRadius = max_dist / 2;  
+}				
+
+template<typename ValType, int Ndims>
+ComputePrecision 
+SPHSmooth<ValType,Ndims>::computeSmoothedValue(const typename SPHTree::coords& c,
+					       computeParticleValue fun)
+{
+  ComputePrecision outputValue = 0;
+  ComputePrecision max_dist = 0;
+  ComputePrecision r3 = smoothRadius * smoothRadius * smoothRadius;
+
+  for (uint32_t i = 0; i < currentNgb; i++)
+    {
+      outputValue += computeWValue(c, *ngb[i], distances[i], fun);
+    }
+
+  return outputValue / r3;
+}
+
+template<typename ValType>
+ComputePrecision interpolateOne(const ValType& t)
+{
+  return 1.0;
+}
+
+// WARNING ! Cell's weight must be 1 !!!
+template<typename ValType, int Ndims>
+ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typename SPHTree::coords& c,
+								    computeParticleValue fun)
+{
+  ComputePrecision outputValue = 0;
+  ComputePrecision max_dist = 0;
+  ComputePrecision weight = 0;
+
+  for (uint32_t i = 0; i < currentNgb; i++)
+    {
+      outputValue += computeWValue(c, *ngb[i], distances[i], fun);
+      weight += computeWValue(c, *ngb[i], distances[i], interpolateOne);
+    }
+
+  return (outputValue == 0) ? 0 : (outputValue / weight);  
+}
+
+template<typename ValType, int Ndims>
+void SPHSmooth<ValType,Ndims>::runForEachNeighbour(runParticleValue fun)
+{
+  for (uint32_t i = 0; i < currentNgb; i++)
+    {
+      fun(ngb[i]);
+    }
+}
+
+
+template<typename ValType, int Ndims>
+void SPHSmooth<ValType,Ndims>::addGridSite(const typename SPHTree::coords& c)
+{
+  ComputePrecision outputValue = 0;
+  ComputePrecision max_dist = 0;
+  
+  ComputePrecision r3 = smoothRadius * smoothRadius * smoothRadius;
+
+  for (uint32_t i = 0; i < currentNgb; i++)
+    {
+      ComputePrecision d = distances[i];
+      SPHCell& cell = *ngb[i];
+      cell.val.weight += getKernel(d/smoothRadius) / r3;
+
+    }
+}
+
+template<typename ValType, int Ndims>
+ComputePrecision 
+SPHSmooth<ValType,Ndims>::getKernel(ComputePrecision x) const
+{
+  // WARNING !!! This is an unnormalized version of the kernel.
+  if (x < 1)
+    return 1 - 1.5 * x * x + 0.75 * x * x * x;
+  else if (x < 2)
+    {
+      ComputePrecision d = 2 - x;
+      return 0.25 * d * d * d;
+    }
+  else
+    return 0;
+}
+
+template<typename ValType, int Ndims>
+bool SPHSmooth<ValType,Ndims>::hasNeighbours() const
+{
+  return (currentNgb != 0);
+}
+
+template<class ValType1, class ValType2, int Ndims>
+bool operator<(const SPHSmooth<ValType1, Ndims>& s1, const SPHSmooth<ValType2, Ndims>& s2)
+{
+  return (s1.getSmoothingLen() < s2.getSmoothingLen());
+}
+
+};

external/cosmotool/src/yorick.hpp

@@ -0,0 +1,196 @@
+#ifndef __YORICK_HELPERS_HPP
+#define __YORICK_HELPERS_HPP
+
+#include "config.hpp"
+#include <stdint.h>
+#include <fstream>
+#include <string>
+
+
+namespace CosmoTool 
+{
+
+  class ProgressiveDoubleOutputImpl
+  {
+  public:
+    virtual ~ProgressiveDoubleOutputImpl();
+    virtual void addDouble(double d) = 0;
+  };
+
+  template<class T>
+  class ProgressiveInputImpl
+  {
+  public:
+    virtual ~ProgressiveInputImpl() {}
+    virtual T read() = 0;
+    virtual void seek(uint32_t *pos) = 0;
+  };
+
+  template<class T>
+  class ProgressiveOutputImpl
+  {
+  public:
+    virtual ~ProgressiveOutputImpl() {}
+    virtual void put(T a) = 0;
+  };
+
+  class ProgressiveDoubleOutput
+  {
+  private:
+    bool initialized;
+    int *ref;
+    ProgressiveDoubleOutputImpl *impl;
+    
+    friend ProgressiveDoubleOutput saveDoubleArrayProgressive(const char *fname, uint32_t *dimList, uint32_t rank);
+    void decRef();
+  public:
+    ProgressiveDoubleOutput();
+    ProgressiveDoubleOutput(ProgressiveDoubleOutputImpl *i);
+    ProgressiveDoubleOutput(const ProgressiveDoubleOutput& o);
+    ~ProgressiveDoubleOutput();
+
+    virtual void addDouble(double a);
+
+    const ProgressiveDoubleOutput& operator=(const ProgressiveDoubleOutput& b); 
+  };
+
+  template<class T>
+  class ProgressiveInput
+  {
+  private:
+    int *ref;
+    ProgressiveInputImpl<T> *impl;
+     
+    void decRef()
+    {
+      if (ref == 0)
+	return;
+      
+      (*ref)--;
+      if (*ref == 0)
+	{
+      	  delete ref;
+	  delete impl;
+	}
+      impl = 0;
+      ref = 0;
+    }
+
+  public:
+    static ProgressiveInput<T> 
+    loadArrayProgressive(const std::string& fname, uint32_t *&dimList, 
+			 uint32_t& rank);
+
+    ProgressiveInput() {
+      impl = 0; 
+      ref = 0; 
+    }
+    ProgressiveInput(ProgressiveInputImpl<T> *i) {
+      impl = i; 
+      ref = new int; 
+      *ref = 1;
+    }
+    ProgressiveInput(const ProgressiveInput<T>& o) {
+      ref = o.ref;
+      impl = o.impl;
+      (*ref)++;
+    }
+    ~ProgressiveInput() {
+      decRef();
+    }
+
+    T read()
+    {
+      return impl->read();
+    }
+
+    void seek(uint32_t *pos)
+    {
+      impl->seek(pos);
+    }
+
+    const ProgressiveInput<T>& operator=(const ProgressiveInput<T>& b)
+    {
+      decRef();
+      ref = b.ref;
+      impl = b.impl;
+      if (ref != 0)
+	(*ref)++;
+      return *this;
+    }
+  };
+
+  template<class T>
+  class ProgressiveOutput
+  {
+  private:
+    int *ref;
+    ProgressiveOutputImpl<T> *impl;
+     
+    void decRef()
+    {
+      if (ref == 0)
+	return;
+      
+      (*ref)--;
+      if (*ref == 0)
+	{
+      	  delete ref;
+	  delete impl;
+	}
+      impl = 0;
+      ref = 0;
+    }
+
+  public:
+    static ProgressiveOutput<T> 
+    saveArrayProgressive(const std::string& fname, uint32_t *dimList, 
+			 uint32_t rank);
+
+    ProgressiveOutput() {
+      impl = 0; 
+      ref = 0; 
+    }
+    ProgressiveOutput(ProgressiveOutputImpl<T> *i) {
+      impl = i; 
+      ref = new int; 
+      *ref = 1;
+    }
+    ProgressiveOutput(const ProgressiveOutput<T>& o) {
+      ref = o.ref;
+      impl = o.impl;
+      (*ref)++;
+    }
+    ~ProgressiveOutput() {
+      decRef();
+    }
+
+    void put(T a)
+    {
+      impl->put(a);
+    }
+
+    const ProgressiveOutput<T>& operator=(const ProgressiveOutput<T>& b)
+    {
+      decRef();
+      ref = b.ref;
+      impl = b.impl;
+      if (ref != 0)
+	(*ref)++;
+      return *this;
+    }
+  };
+
+  template<typename T>
+  void saveArray(const std::string& fname,
+		 T *array, uint32_t *dimList, uint32_t rank);
+
+  template<typename T>
+  void loadArray(const std::string& fname,
+		 T*& array, uint32_t *& dimList, uint32_t& rank)
+    throw (NoSuchFileException);
+
+  ProgressiveDoubleOutput saveDoubleArrayProgressive(const char *fname, uint32_t *dimList, uint32_t rank);  
+};
+
+#endif

external/cosmotool/src/yorick_nc3.cpp

@@ -0,0 +1,274 @@
+#include "ctool_netcdf_ver.hpp"
+#include "config.hpp"
+#ifdef NETCDFCPP4
+#include <netcdf>
+using namespace netCDF
+#else
+#include <netcdfcpp.h>
+#endif
+#include <fstream>
+#include "yorick.hpp"
+#include <assert.h>
+
+using namespace CosmoTool;
+using namespace std;
+
+class NetCDF_handle
+{
+public:
+  NcFile *outFile;
+  NcVar *curVar;
+  long *curPos;
+  long *counts;
+  long *dimList;
+  uint32_t rank;
+  
+  NetCDF_handle(NcFile *f, NcVar *v, long *dimList, uint32_t rank);
+  virtual ~NetCDF_handle();
+};
+
+NetCDF_handle::NetCDF_handle(NcFile *f, NcVar *v, long *dimList, uint32_t rank)
+{
+  this->outFile = f;
+  this->curVar = v;
+  this->dimList = dimList;
+  this->rank = rank;
+  this->counts = new long[rank];
+  this->curPos = new long[rank];
+
+  for (long i = 0; i < rank; i++)
+    this->curPos[i] = 0;
+
+  for (long i = 0; i < rank; i++)
+    this->counts[i] = 1;
+}
+
+NetCDF_handle::~NetCDF_handle()
+{
+  delete[] dimList;
+  delete outFile;
+}
+
+template<typename T>
+class InputGenCDF: public NetCDF_handle, public ProgressiveInputImpl<T>
+{
+public:
+  InputGenCDF(NcFile *f, NcVar *v, long *dimList, uint32_t rank)
+    : NetCDF_handle(f,v,dimList,rank)
+  {}
+  virtual ~InputGenCDF() {} 
+
+  virtual T read()
+  {
+    T a;
+
+    curVar->set_cur(curPos);
+    curVar->get(&a, counts);
+    
+    curPos[rank-1]++;
+    for (long i = rank-1; i >= 1; i--)
+      {
+	if (curPos[i] == dimList[i])
+	  {
+	    curPos[i-1]++;
+	  curPos[i] = 0;
+	}
+      }    
+    return a;
+  }
+
+  virtual void seek(uint32_t *pos)
+  {
+    for (long i = rank-1; i >= 0; i--)
+      curPos[i] = pos[rank-1-i];
+  }
+};
+
+template<typename T>
+class OutputGenCDF: public NetCDF_handle, public ProgressiveOutputImpl<T>
+{
+public:
+  OutputGenCDF(NcFile *f, NcVar *v, long *dimList, uint32_t rank)
+    : NetCDF_handle(f,v,dimList,rank)
+  {}
+  virtual ~OutputGenCDF() {} 
+
+  virtual void put(T a)
+  {
+    curVar->set_cur(curPos);
+    curVar->put(&a, counts);
+    
+    curPos[rank-1]++;
+    for (long i = rank-1; i >= 1; i--)
+      {
+	if (curPos[i] == dimList[i])
+	  {
+	    curPos[i-1]++;
+	  curPos[i] = 0;
+	}
+      }    
+  }
+};
+
+template<typename T>
+class NetCDF_type
+{
+public:
+  static const NcType t = (NcType)-1;
+};
+
+template<>
+class NetCDF_type<int>
+{
+public:
+  static const NcType t = ncInt;
+};
+
+template<>
+class NetCDF_type<unsigned int>
+{
+public:
+  static const NcType t = ncInt;
+};
+
+template<>
+class NetCDF_type<float>
+{
+public:
+  static const NcType t = ncFloat;
+};
+
+template<>
+class NetCDF_type<double>
+{
+public:
+  static const NcType t = ncDouble;
+};
+
+namespace CosmoTool {
+  template<typename T>
+  ProgressiveOutput<T>
+  ProgressiveOutput<T>::saveArrayProgressive(const std::string& fname, uint32_t *dimList,
+					     uint32_t rank)
+  {
+    NcFile *f = new NcFile(fname.c_str(), NcFile::Replace);
+    
+    assert(f->is_valid());
+    
+    const NcDim **dimArray = new const NcDim *[rank];
+    for (uint32_t i = 0; i < rank; i++)
+      {
+	char dimName[255];
+	
+	sprintf(dimName, "dim%d", i);
+	dimArray[i] = f->add_dim(dimName, dimList[rank-1-i]);
+      }
+    
+    NcVar *v = f->add_var("array", NetCDF_type<T>::t, rank, dimArray);  
+
+    long *ldimList = new long[rank];
+
+    for (uint32_t i = 0; i < rank; i++)
+      ldimList[rank-1-i] = dimList[i];
+
+    OutputGenCDF<T> *impl = new OutputGenCDF<T>(f, v, ldimList, rank);
+    return ProgressiveOutput<T>(impl);   
+  }
+
+  template<typename T>
+  ProgressiveInput<T>
+  ProgressiveInput<T>::loadArrayProgressive(const std::string& fname, uint32_t *&dimList,
+					      uint32_t& rank)
+  {
+    NcFile *f = new NcFile(fname.c_str(), NcFile::ReadOnly);
+
+    assert(f->is_valid());
+
+    NcVar *v = f->get_var("array");
+
+    rank = v->num_dims();
+    long *ldimList = v->edges();
+    dimList = new uint32_t[rank];
+    for (uint32_t i = 0; i < rank; i++)
+      {
+	dimList[rank-i-1] = ldimList[i];
+      }
+    InputGenCDF<T> *impl = new InputGenCDF<T>(f, v, ldimList, rank);
+
+    return ProgressiveInput<T>(impl);
+  }
+
+  template<typename T>
+  void saveArray(const std::string& fname,
+		 T *array, uint32_t *dimList, uint32_t rank)
+  {
+    NcFile f(fname.c_str(), NcFile::Replace);
+    
+    assert(f.is_valid());
+    
+    const NcDim **dimArray = new const NcDim *[rank];
+    for (uint32_t i = 0; i < rank; i++)
+      {
+	char dimName[255];
+	
+	sprintf(dimName, "dim%d", i);
+	dimArray[i] = f.add_dim(dimName, dimList[i]);
+      }
+
+    NcVar *v = f.add_var("array", NetCDF_type<T>::t, rank, dimArray);
+    
+    long *edge = v->edges();
+    v->put(array, edge);
+    delete[] edge;  
+  }
+
+  template<typename T>
+  void loadArray(const std::string& fname,
+		 T*&array, uint32_t *&dimList, uint32_t& rank)
+	throw (NoSuchFileException)
+  {
+    NcFile f(fname.c_str(), NcFile::ReadOnly);
+    
+    if (!f.is_valid())
+      throw NoSuchFileException(fname);
+    
+    NcVar *v = f.get_var("array"); 
+    rank = v->num_dims();
+    long *edge = v->edges();
+    uint32_t fullSize = 1;
+    dimList = new uint32_t[rank];
+    for (int i = 0; i < rank; i++)
+      {
+	dimList[i] = edge[i];
+	fullSize *= edge[i];
+      }
+    if (fullSize != 0) {
+      array = new T[fullSize];
+      v->get(array, edge);
+    }
+    delete[] edge;   
+  }
+
+  template class ProgressiveInput<int>;
+  template class ProgressiveInput<float>;
+  template class ProgressiveInput<double>;
+
+  template class ProgressiveOutput<int>;
+  template class ProgressiveOutput<float>;
+  template class ProgressiveOutput<double>;
+
+  template void loadArray<int>(const std::string& fname,
+			       int*& array, uint32_t *&dimList, uint32_t& rank);
+  template void loadArray<float>(const std::string& fname,
+				 float*& array, uint32_t *&dimList, uint32_t& rank);
+  template void loadArray<double>(const std::string& fname,
+				  double*& array, uint32_t *&dimList, uint32_t& rank);
+
+  template void saveArray<int>(const std::string& fname,
+			       int *array, uint32_t *dimList, uint32_t rank);
+  template void saveArray<float>(const std::string& fname,
+				 float *array, uint32_t *dimList, uint32_t rank);
+  template void saveArray<double>(const std::string& fname,
+				  double *array, uint32_t *dimList, uint32_t rank);
+  
+}

external/cosmotool/src/yorick_nc4.cpp

@@ -0,0 +1,237 @@
+#include "ctool_netcdf_ver.hpp"
+#include "config.hpp"
+#include <netcdf>
+using namespace netCDF;
+#include <fstream>
+#include "yorick.hpp"
+#include <assert.h>
+
+using namespace CosmoTool;
+using namespace std;
+
+class NetCDF_handle
+{
+public:
+  NcFile *outFile;
+  NcVar curVar;
+  vector<size_t> curPos;
+  vector<size_t> counts;
+  vector<NcDim> dimList;
+  uint32_t rank;
+  
+  NetCDF_handle(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank);
+  virtual ~NetCDF_handle();
+};
+
+NetCDF_handle::NetCDF_handle(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank)
+{
+  this->outFile = f;
+  this->curVar = v;
+  this->dimList = dimList;
+  this->rank = rank;
+  this->counts.resize(rank);
+  this->curPos.resize(rank);
+
+  for (long i = 0; i < rank; i++)
+    this->curPos[i] = 0;
+
+  for (long i = 0; i < rank; i++)
+    this->counts[i] = 1;
+}
+
+NetCDF_handle::~NetCDF_handle()
+{
+  delete outFile;
+}
+
+template<typename T>
+class InputGenCDF: public NetCDF_handle, public ProgressiveInputImpl<T>
+{
+public:
+  InputGenCDF(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank)
+    : NetCDF_handle(f,v,dimList,rank)
+  {}
+  virtual ~InputGenCDF() {} 
+
+  virtual T read()
+  {
+    T a;
+
+    curVar.getVar(curPos, counts, &a);
+    
+    curPos[rank-1]++;
+    for (long i = rank-1; i >= 1; i--)
+      {
+	if (curPos[i] == dimList[i].getSize())
+	  {
+	    curPos[i-1]++;
+	  curPos[i] = 0;
+	}
+      }    
+    return a;
+  }
+
+  virtual void seek(uint32_t *pos)
+  {
+    for (long i = rank-1; i >= 0; i--)
+      curPos[i] = pos[rank-1-i];
+  }
+};
+
+template<typename T>
+class OutputGenCDF: public NetCDF_handle, public ProgressiveOutputImpl<T>
+{
+public:
+  OutputGenCDF(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank)
+    : NetCDF_handle(f,v,dimList,rank)
+  {}
+  virtual ~OutputGenCDF() {} 
+
+  virtual void put(T a)
+  {
+    curVar.putVar(curPos, counts, &a);
+    
+    curPos[rank-1]++;
+    for (long i = rank-1; i >= 1; i--)
+      {
+	if (curPos[i] == dimList[i].getSize())
+	  {
+	    curPos[i-1]++;
+	  curPos[i] = 0;
+	}
+      }    
+  }
+};
+
+template<typename T> NcType& get_NetCDF_type();
+
+#define IMPL_TYPE(T,ncT) \
+template<> \
+NcType& get_NetCDF_type<T>() \
+{ \
+  return ncT; \
+}
+
+IMPL_TYPE(int,ncInt);
+IMPL_TYPE(unsigned int,ncInt);
+IMPL_TYPE(double,ncDouble);
+IMPL_TYPE(float,ncFloat);
+
+namespace CosmoTool {
+  template<typename T>
+  ProgressiveOutput<T>
+  ProgressiveOutput<T>::saveArrayProgressive(const std::string& fname, uint32_t *dimList,
+					     uint32_t rank)
+  {
+    NcFile *f = new NcFile(fname, NcFile::replace);
+    
+    vector<NcDim> dimArray;
+    for (uint32_t i = 0; i < rank; i++)
+      {
+	char dimName[255];
+	
+	sprintf(dimName, "dim%d", i);
+	dimArray.push_back(f->addDim(dimName, dimList[rank-1-i]));
+      }
+    
+    NcVar v = f->addVar("array", get_NetCDF_type<T>(), dimArray);  
+
+    vector<NcDim> ldimList;
+
+    for (uint32_t i = 0; i < rank; i++)
+      ldimList.push_back(dimArray[rank-1-i]);
+
+    OutputGenCDF<T> *impl = new OutputGenCDF<T>(f, v, ldimList, rank);
+    return ProgressiveOutput<T>(impl);   
+  }
+
+  template<typename T>
+  ProgressiveInput<T>
+  ProgressiveInput<T>::loadArrayProgressive(const std::string& fname, uint32_t *&dimList,
+					      uint32_t& rank)
+  {
+    NcFile *f = new NcFile(fname, NcFile::read);
+
+    NcVar v = f->getVar("array");
+
+    rank = v.getDimCount();
+    vector<NcDim> vdimlist = v.getDims();
+    dimList = new uint32_t[rank];
+    for (uint32_t i = 0; i < rank; i++)
+      {
+	dimList[rank-i-1] = vdimlist[i].getSize();
+      }
+    InputGenCDF<T> *impl = new InputGenCDF<T>(f, v, vdimlist, rank);
+
+    return ProgressiveInput<T>(impl);
+  }
+
+  template<typename T>
+  void saveArray(const std::string& fname,
+		 T *array, uint32_t *dimList, uint32_t rank)
+  {
+    NcFile f(fname.c_str(), NcFile::replace);
+    
+    vector<NcDim> dimArray;
+    for (uint32_t i = 0; i < rank; i++)
+      {
+	char dimName[255];
+	
+	sprintf(dimName, "dim%d", i);
+	dimArray.push_back(f.addDim(dimName, dimList[i]));
+      }
+
+    NcVar v = f.addVar("array", get_NetCDF_type<T>(), dimArray);
+    
+    v.putVar(array);
+  }
+
+  template<typename T>
+  void loadArray(const std::string& fname,
+		 T*&array, uint32_t *&dimList, uint32_t& rank)
+	throw (NoSuchFileException)
+  {
+    NcFile f(fname.c_str(), NcFile::read);
+    
+    //if (!f.is_valid())
+    //  throw NoSuchFileException(fname);
+    
+    NcVar v = f.getVar("array"); 
+    vector<NcDim> dims = v.getDims();
+    rank = v.getDimCount();
+    uint32_t fullSize = 1;
+    dimList = new uint32_t[rank];
+    for (int i = 0; i < rank; i++)
+      {
+	dimList[i] = dims[i].getSize();
+	fullSize *= dimList[i];
+      }
+    if (fullSize != 0) {
+      array = new T[fullSize];
+      v.getVar(array);
+    }
+  }
+
+  template class ProgressiveInput<int>;
+  template class ProgressiveInput<float>;
+  template class ProgressiveInput<double>;
+
+  template class ProgressiveOutput<int>;
+  template class ProgressiveOutput<float>;
+  template class ProgressiveOutput<double>;
+
+  template void loadArray<int>(const std::string& fname,
+			       int*& array, uint32_t *&dimList, uint32_t& rank);
+  template void loadArray<float>(const std::string& fname,
+				 float*& array, uint32_t *&dimList, uint32_t& rank);
+  template void loadArray<double>(const std::string& fname,
+				  double*& array, uint32_t *&dimList, uint32_t& rank);
+
+  template void saveArray<int>(const std::string& fname,
+			       int *array, uint32_t *dimList, uint32_t rank);
+  template void saveArray<float>(const std::string& fname,
+				 float *array, uint32_t *dimList, uint32_t rank);
+  template void saveArray<double>(const std::string& fname,
+				  double *array, uint32_t *dimList, uint32_t rank);
+  
+}