Merge branch 'master' of file:///home/lavaux/Dropbox/gitRoot/CosmoToolbox

src/CMakeLists.txt

@@ -6,11 +6,19 @@ SET(CosmoTool_SRCS
 	loadRamses.cpp
 	octTree.cpp
 	powerSpectrum.cpp
-	yorick.cpp
 	miniargs.cpp
 	growthFactor.cpp
 )
 
+IF(FOUND_NETCDF3)
+	SET(CosmoTool_SRCS ${CosmoTool_SRCS} yorick_nc3.cpp)
+ELSE(FOUND_NETCDF3)
+IF(FOUND_NETCDF4)
+	SET(CosmoTool_SRCS ${CosmoTool_SRCS} yorick_nc4.cpp)
+ENDIF(FOUND_NETCDF4)
+ENDIF(FOUND_NETCDF3)
+
+
 if (HDF5_FOUND)
   set(CosmoTool_SRCS ${CosmoTool_SRCS}
     h5_readFlash.cpp
@@ -46,14 +54,16 @@ SET(CosmoTool_SRCS ${CosmoTool_SRCS}
 	growthFactor.hpp
 )
 
-include_directories(${GSL_INCLUDE_PATH} ${NETCDF_INCLUDE_PATH})
+include_directories(${GSL_INCLUDE_PATH} ${NETCDF_INCLUDE_PATH} ${NETCDFCPP_INCLUDE_PATH} ${CMAKE_BINARY_DIR}/src)
 
-set(CosmoTool_LIBS  ${NETCDF_LIBRARY} ${NETCDFCPP_LIBRARY} ${GSL_LIBRARIES})
+set(CosmoTool_LIBS  ${NETCDFCPP_LIBRARY} ${NETCDF_LIBRARY} ${GSL_LIBRARIES})
 if (HDF5_FOUND)
   set(CosmoTool_LIBS ${CosmoTool_LIBS} ${HDF5_CXX_LIBRARIES} ${HDF5_LIBRARIES})
   include_directories(${HDF5_INCLUDE_DIRS})
 endif (HDF5_FOUND)
 
+set(CosmoTool_LIBS ${CosmoTool_LIBS} PARENT_SCOPE)
+
 if (BUILD_SHARED_LIBS)
   add_library(CosmoTool SHARED ${CosmoTool_SRCS})
   target_link_libraries(CosmoTool ${CosmoTool_LIBS})
@@ -62,7 +72,7 @@ if (BUILD_SHARED_LIBS)
      add_library(CosmoTool_static STATIC ${CosmoTool_SRCS})
   endif(BUILD_STATIC_LIBS)
 else (BUILD_SHARED_LIBS)
-  add_library(CosmoTool STATIC ${CosmoTool_LIBS})
+  add_library(CosmoTool STATIC ${CosmoTool_SRCS})
 endif (BUILD_SHARED_LIBS)
 
 install(TARGETS CosmoTool 

src/bsp_simple.hpp

@@ -154,6 +154,8 @@ namespace CosmoTool
 	f.data = data;
 	insert(f);
       }
+
+      bool inside(const typename space_t::coord_t& p) const;
     };
 
   };

src/bsp_simple.tcc

@@ -98,7 +98,19 @@ namespace CosmoTool
 	*(*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
+    }
 
   };
 

src/fourier/base_types.hpp

@@ -0,0 +1,136 @@
+#ifndef __BASE_FOURIER_TYPES_HPP
+#define __BASE_FOURIER_TYPES_HPP
+
+#include <gsl/gsl_rng.h>
+#include <boost/shared_ptr.hpp>
+#include <string>
+#include <Eigen/Dense>
+#include <complex>
+#include <exception>
+
+namespace CosmoTool
+{
+  class IncompatibleMap: virtual std::exception {};
+
+  template<typename T> class FourierMap;
+
+  template<typename T>
+  class SpectrumFunction
+  {
+  protected:
+    SpectrumFunction() {}
+  public:
+    typedef Eigen::Array<T, 1, Eigen::Dynamic> VecType;
+    typedef Eigen::Map<VecType, Eigen::Aligned> MapType;
+    typedef Eigen::Map<VecType const, Eigen::Aligned> ConstMapType;
+    typedef FourierMap<std::complex<T> > FourierMapType;
+
+    virtual boost::shared_ptr<FourierMapType> 
+       newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0;
+
+    virtual void mul(FourierMap<std::complex<T> >& m) const = 0;
+  };
+
+
+  template<typename T>
+  class FourierMap
+  {
+  protected:
+    FourierMap() {}
+    
+  public:
+    typedef Eigen::Array<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 sqrt()
+    {
+      MapType m = eigen();
+      m = m.sqrt();
+    }
+ 
+    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 T dot_product(const FourierMap<T>& second) const
+       throw(std::bad_cast) = 0;
+
+    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

src/fourier/euclidian.hpp

@@ -0,0 +1,407 @@
+#ifndef __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
+#define __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
+
+#include <boost/function.hpp>
+#include <vector>
+#include <boost/shared_ptr.hpp>
+#include <gsl/gsl_randist.h>
+#include "base_types.hpp"
+#include "fft/fftw_calls.hpp"
+#include "../algo.hpp"
+
+namespace CosmoTool
+{
+  template<typename T>
+  class EuclidianSpectrum_1D: public SpectrumFunction<T>
+  {
+  public:
+    typedef boost::function1<T, T> Function;
+  protected:
+    Function f;
+  public:
+    typedef typename SpectrumFunction<T>::FourierMapType FourierMapType;
+    typedef boost::shared_ptr<FourierMapType> ptr_map;
+
+    EuclidianSpectrum_1D(Function P)
+      : f(P)
+    {
+    }
+
+    ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const;
+
+    void mul(FourierMap<std::complex<T> >& m) const;
+  };
+  
+  template<typename T>
+  class EuclidianFourierMapBase: public FourierMap<T>
+  {
+  public:
+    typedef std::vector<int> DimArray;
+  private:
+    boost::shared_ptr<T> m_data;
+    DimArray m_dims;
+    long m_size;
+  public:
+    
+    EuclidianFourierMapBase(boost::shared_ptr<T> indata, const DimArray& 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 ~EuclidianFourierMapBase()
+    {
+    }
+
+    const DimArray& getDims() const { return m_dims; }
+
+    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;
+    }
+
+  };
+
+  template<typename T>
+  class EuclidianFourierMapReal: public EuclidianFourierMapBase<T>
+  {
+  public:
+    typedef typename EuclidianFourierMapBase<T>::DimArray DimArray;
+
+    EuclidianFourierMapReal(boost::shared_ptr<T> indata, const DimArray& indims)
+      : EuclidianFourierMapBase<T>(indata, indims) 
+    {}
+
+    virtual FourierMap<T> *mimick() const
+    {
+      return new EuclidianFourierMapReal<T>(
+		      boost::shared_ptr<T>((T *)fftw_malloc(sizeof(T)*this->size()), 
+					   std::ptr_fun(fftw_free)),
+		      this->getDims());
+    }
+
+    virtual T dot_product(const FourierMap<T>& other) const
+      throw(std::bad_cast)
+    {
+      const EuclidianFourierMapReal<T>& m2 = dynamic_cast<const EuclidianFourierMapReal<T>&>(other);
+      if (this->size() != m2.size())
+        throw std::bad_cast();
+
+      return (this->eigen()*m2.eigen()).sum();
+    }
+  };
+
+  template<typename T>
+  class EuclidianFourierMapComplex: public EuclidianFourierMapBase<std::complex<T> >
+  {
+  protected:
+    typedef boost::shared_ptr<std::complex<T> > ptr_t;
+    std::vector<double> delta_k;
+    long plane_size;
+  public:
+    typedef typename EuclidianFourierMapBase<std::complex<T> >::DimArray DimArray;
+
+    EuclidianFourierMapComplex(ptr_t indata,
+			       const DimArray& indims, 
+			       const std::vector<double>& dk)
+      : EuclidianFourierMapBase<std::complex<T> >(indata, indims), delta_k(dk)
+    {
+      assert(dk.size() == indims.size()); 
+      plane_size = 1;
+      for (int q = 1; q < indims.size(); q++)
+	plane_size *= indims[q];
+    }
+
+    virtual FourierMap<std::complex<T> > *mimick() const
+    {
+      return 
+	new EuclidianFourierMapComplex<T>(
+		   ptr_t((std::complex<T> *)
+			 fftw_malloc(sizeof(std::complex<T>)*this->size()), 
+			 std::ptr_fun(fftw_free)),
+		   this->getDims(),
+		   this->delta_k);
+    }
+
+    template<typename Array>
+    double get_K(const Array& ik)
+    {
+      const DimArray& dims = this->getDims();
+      assert(ik.size() == dims.size());
+      double k2 = 0;
+      for (int q = 0; q < ik.size(); q++)
+	{
+	  int dk = ik;
+
+	  if (dk > dims[q]/2)
+	    dk = dk - dims[q];
+	  
+	  k2 += CosmoTool::square(delta_k[q]*dk);
+	}
+      return std::sqrt(k2);
+    }
+
+    double get_K(long p)
+    {
+      const DimArray& dims = this->getDims();
+      DimArray d(delta_k.size());
+      for (int q = 0; q < d.size(); q++)
+       {
+         d[q] = p%dims[q];
+         p = (p-d[q])/dims[q];
+       }
+      return get_K(d);
+    }
+
+    virtual std::complex<T> dot_product(const FourierMap<std::complex<T> >& other) const
+      throw(std::bad_cast)
+    {
+      const EuclidianFourierMapComplex<T>& m2 = dynamic_cast<const EuclidianFourierMapComplex<T>&>(other);
+      if (this->size() != m2.size())
+        throw std::bad_cast();
+      
+      const std::complex<T> *d1 = this->data();
+      const std::complex<T> *d2 = m2.data();
+      const DimArray& dims = this->getDims();
+      int N0 = dims[0];
+      std::complex<T> result = 0;
+
+      for (long q0 = 1; q0 < N0-1; q0++)
+	{
+	  for (long p = 0; p < plane_size; p++)
+	    {
+	      long idx = q0+N0*p;
+	      assert(idx < this->size());
+	      result += 2*(conj(d1[idx]) * d2[idx]).real();
+	    }
+	}
+      for (long p = 0; p < plane_size; p++)
+	{
+	  long q0 = N0*p, q1 = (p+1)*N0-1;
+	  result += conj(d1[q0]) * d2[q0];
+	  result += conj(d1[q1]) * d2[q1];
+	}
+      return result;
+    }
+
+  };
+  
+  template<typename T>
+  class EuclidianFourierTransform: public FourierTransform<T>
+  {
+  public:
+    typedef typename EuclidianFourierMapBase<T>::DimArray DimArray;
+  private:
+    typedef FFTW_Calls<T> calls;
+    EuclidianFourierMapReal<T> *realMap;
+    EuclidianFourierMapComplex<T> *fourierMap;
+    typename calls::plan_type m_analysis, m_synthesis;
+    double volume;
+    long N;
+    DimArray m_dims, m_dims_hc;
+    std::vector<double> m_L;
+  public:
+    EuclidianFourierTransform(const DimArray& dims, const std::vector<double>& L)
+    {
+      assert(L.size() == dims.size());
+      std::vector<double> dk(L.size());
+      
+      m_dims = dims;
+      m_dims_hc = dims;
+      m_dims_hc[0] = dims[0]/2+1;
+      m_L = L;
+
+      N = 1;
+      volume = 1;
+      for (int i = 0; i < dims.size(); i++)
+	{
+	  N *= dims[i];
+	  volume *= L[i];
+	  dk[i] = 2*M_PI/L[i];
+	}
+
+      realMap = new EuclidianFourierMapReal<T>(
+		   boost::shared_ptr<T>(calls::alloc_real(N),
+					std::ptr_fun(calls::free)),
+		   m_dims);
+      fourierMap = new EuclidianFourierMapComplex<T>(
+		   boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(N),
+						       std::ptr_fun(calls::free)), 
+		   m_dims_hc, dk);
+      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() {}
+  };
+
+
+  template<typename T>
+  typename EuclidianSpectrum_1D<T>::ptr_map
+   EuclidianSpectrum_1D<T>::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const
+  {
+    typedef EuclidianFourierMapComplex<T> MapT;
+    typedef typename MapT::DimArray DimArray;
+
+    MapT& m_c = dynamic_cast<MapT&>(like_map);
+    MapT *rand_map = m_c.mimick();
+    std::complex<T> *d = rand_map->data();
+    long idx;
+    const DimArray& dims = rand_map->getDims();
+    long plane_size;
+
+    for (long p = 1; p < m_c.size(); p++)
+      {
+	double A_k = std::sqrt(0.5*f(rand_map->get_K(p)));
+	d[p] = std::complex<T>(gsl_ran_gaussian(rng, A_k),
+			       gsl_ran_gaussian(rng, A_k));
+      }
+    // Generate the mean value
+    d[0] = std::complex<T>(std::sqrt(f(0)), 0);
+    // Correct the Nyquist plane
+    idx = dims[0]-1; // Stick to the last element of the first dimension
+    // 1D is special case
+    if (dims.size() == 1)
+      {
+	d[idx] = std::complex<T>(d[idx].real() + d[idx].imag(), 0);
+	return boost::shared_ptr<EuclidianSpectrum_1D<T>::FourierMapType>(rand_map);
+      }
+
+    plane_size = 1;
+    for (int q = 1; q < dims.size(); q++)
+      {
+	plane_size *= dims[q];
+      }
+
+    for (long p = 1; p < plane_size/2; p++)
+      {
+	long q = (p+1)*dims[0]-1;
+	long q2 = (plane_size-p+1)*dims[0]-1;
+	assert(q < plane_size*dims[0]);
+	assert(q2 < plane_size*dims[0]);
+	d[q] = conj(d[q2]);
+      }
+    long q = dims[0];
+    d[q] = std::complex<T>(d[q].real() + d[q].imag());
+  }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul(FourierMap<std::complex<T> >& m) const
+  {
+    EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m);
+    std::complex<T> *d = m.data();
+
+    for (long p = 0; p < m_c.size(); p++)
+      d[p] *= f(m.get_K(p));
+  }
+};
+
+#endif

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

src/fourier/healpix.hpp

@@ -0,0 +1,268 @@
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_HPP
+#define __COSMOTOOL_FOURIER_HEALPIX_HPP
+
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <cmath>
+#include <vector>
+#include <boost/bind.hpp>
+#include <boost/shared_ptr.hpp>
+#include <exception>
+#include "base_types.hpp"
+#include <sharp_lowlevel.h>
+#include <sharp_geomhelpers.h>
+#include <sharp_almhelpers.h>
+
+namespace CosmoTool
+{
+  template<typename T>
+  class HealpixSpectrum: public SpectrumFunction<T>
+  {
+  protected:
+    std::vector<T> cls;
+  public:
+    typedef typename SpectrumFunction<T>::FourierMapType FourierMapType;
+    typedef boost::shared_ptr<FourierMapType> ptr_map;
+
+
+    HealpixSpectrum(long Lmax)
+     : cls(Lmax+1) {}
+
+    T *data() { return &cls[0]; }
+    const T *data() const { return &cls[0]; }
+    long size() const { return cls.size(); }
+
+    ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0;
+  };    
+
+  template<typename T>
+  class HealpixFourierMap: public FourierMap<T>
+  {
+  private:
+    T *m_data;
+    long Npix, m_Nside;
+    Eigen::aligned_allocator<T> alloc;
+  public:
+    HealpixFourierMap(long nSide)
+      : Npix(12*nSide*nSide), m_Nside(nSide)
+    {
+      m_data = alloc.allocate(Npix);
+    }
+
+    virtual ~HealpixFourierMap()
+    {
+      alloc.deallocate(m_data, Npix);
+    }
+
+    long Nside() const { return m_Nside; }
+    virtual const T* data() const { return m_data; }
+    virtual T *data() { return m_data; }
+    virtual long size() const { return Npix; }
+
+    virtual T dot_product(const FourierMap<T>& other) const
+      throw(std::bad_cast)
+    {
+      typedef typename FourierMap<T>::MapType MapType;
+
+      const HealpixFourierMap<T>& mfm = dynamic_cast<const HealpixFourierMap<T>&>(other);
+      if (Npix != mfm.size())
+        throw std::bad_cast();
+
+      MapType m1(m_data, Npix);
+      MapType m2(mfm.m_data, mfm.Npix);
+      
+      return (m1*m2).sum();
+    }
+    
+    virtual FourierMap<T> *mimick() const
+    {
+      return new HealpixFourierMap<T>(m_Nside);
+    }
+  };
+
+
+  template<typename T>
+  class HealpixFourierALM: public FourierMap<std::complex<T> >
+  {
+  private:
+    std::complex<T> *alms;
+    long m_size;
+    long Lmax_, Mmax_, TVal_;
+    Eigen::aligned_allocator<std::complex<T> > alloc;
+  public:
+    typedef unsigned long LType;
+
+    LType Lmax() const { return Lmax_; }
+    LType Mmax() const { return Mmax_; }
+
+    LType Num_Alms() const
+    {
+      return ((Mmax_+1)*(Mmax_+2))/2 + (Mmax_+1)*(Lmax_-Mmax_);
+    }
+
+    LType index_l0(LType m) const
+    {
+      return ((m*(TVal_-m))/2);
+    }
+
+    LType index(LType l, LType m) const
+    {
+      return index_l0(m) + l;
+    }
+    
+    HealpixFourierALM(LType lmax, LType mmax)
+      : Lmax_(lmax), Mmax_(mmax), TVal_(2*lmax+1)
+    {
+      m_size = Num_Alms();
+      alms = alloc.allocate(m_size);
+    }
+
+    virtual ~HealpixFourierALM()
+    {
+      alloc.deallocate(alms, m_size);
+    }
+
+    virtual const std::complex<T>* data() const { return alms; }
+    virtual std::complex<T> * data() { return alms;} 
+    virtual long size() const { return m_size; }
+
+    virtual FourierMap<std::complex<T> > *mimick() const
+    {
+      return new HealpixFourierALM<T>(Lmax_, Mmax_);
+    }
+
+    virtual std::complex<T> dot_product(const FourierMap<std::complex<T> >& other) const
+      throw(std::bad_cast)
+    {
+      const HealpixFourierALM<T>& mfm = dynamic_cast<const HealpixFourierALM<T>&>(other);
+      typedef typename FourierMap<std::complex<T> >::MapType MapType;
+      std::complex<T> S;
+
+      if (m_size != mfm.m_size)
+        throw std::bad_cast();
+
+      MapType m1(alms, m_size);
+      MapType m2(mfm.alms, mfm.m_size);
+      
+      S = (m1.block(0,0,1,Lmax_+1).conjugate() * m2.block(0,0,1,Lmax_+1)).sum();
+      S += std::complex<T>(2,0)*(m1.block(0,1+Lmax_,1,m_size-1-Lmax_).conjugate() * m2.block(0,1+Lmax_,1,m_size-1-Lmax_)).sum();
+      return S;
+    }
+  };
+
+  template<typename T> struct HealpixJobHelper__ {};
+
+  template<> struct HealpixJobHelper__<double>
+    { enum {val=1}; };
+
+  template<> struct HealpixJobHelper__<float>
+    { enum {val=0}; };
+
+
+  template<typename T>
+  class HealpixFourierTransform: public FourierTransform<T>
+  {
+  private:
+    sharp_alm_info *ainfo;
+    sharp_geom_info *ginfo;
+    HealpixFourierMap<T> realMap;
+    HealpixFourierALM<T> fourierMap;
+    int m_iterate;
+  public:
+    HealpixFourierTransform(long nSide, long Lmax, long Mmax, int iterate = 0)
+      : realMap(nSide), fourierMap(Lmax, Mmax), ainfo(0), ginfo(0), m_iterate(iterate)
+    {
+      sharp_make_healpix_geom_info (nSide, 1, &ginfo);
+      sharp_make_triangular_alm_info (Lmax, Mmax, 1, &ainfo);
+    }
+
+    virtual ~HealpixFourierTransform()
+    {
+      sharp_destroy_geom_info(ginfo);
+      sharp_destroy_alm_info(ainfo);
+    }
+
+    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()
+    {
+      void *aptr=reinterpret_cast<void *>(fourierMap.data()), *mptr=reinterpret_cast<void *>(realMap.data());
+
+      sharp_execute (SHARP_MAP2ALM, 0, 0, &aptr, &mptr, ginfo, ainfo, 1,
+        HealpixJobHelper__<T>::val,0,0,0);
+      for (int i = 0; i < m_iterate; i++)
+        {
+          HealpixFourierMap<T> tmp_map(realMap.Nside());          
+          void *tmp_ptr=reinterpret_cast<void *>(tmp_map.data());
+          typename HealpixFourierMap<T>::MapType m0 = tmp_map.eigen();
+          typename HealpixFourierMap<T>::MapType m1 = realMap.eigen();
+
+          sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &tmp_ptr, ginfo, ainfo, 1,
+            HealpixJobHelper__<T>::val,0,0,0);
+          m0 = m1 - m0;
+          sharp_execute (SHARP_MAP2ALM, 0, 1, &aptr, &tmp_ptr, ginfo, ainfo, 1,
+            HealpixJobHelper__<T>::val,0,0,0);
+        }
+    }
+
+    virtual void synthesis()
+    {
+      void *aptr=reinterpret_cast<void *>(fourierMap.data()), *mptr=reinterpret_cast<void *>(realMap.data());
+
+      sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &mptr, ginfo, ainfo, 1,
+        HealpixJobHelper__<T>::val,0,0,0);
+    }
+
+    virtual void analysis_conjugate()
+    {
+      synthesis();
+      realMap.scale(4*M_PI/realMap.size());
+    }
+
+    virtual void synthesis_conjugate()
+    {
+      analysis();
+      fourierMap.scale(realMap.size()/(4*M_PI));
+    }
+    
+  };
+
+  template<typename T>
+  typename HealpixSpectrum<T>::ptr_map 
+    HealpixSpectrum<T>::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const
+  { 
+    const HealpixFourierALM<T>& alms = dynamic_cast<const HealpixFourierALM<T>&>(like_map);
+    HealpixFourierALM<T> *new_alms;
+    ptr_map r(new_alms = new HealpixFourierALM<T>(alms.Lmax(), alms.Mmax()));
+    long lmaxGen = std::min(cls.size()-1, alms.Lmax());
+    std::complex<T> *new_data = new_alms->data();
+
+    for (long l = 0; l < lmaxGen; l++)
+      {
+        double Al = std::sqrt(cls[l]);
+
+        new_data[alms.index(l,0)] = gsl_ran_gaussian(rng, Al);
+        Al *= M_SQRT1_2;
+        for (long m = 1; m < alms.Mmax(); m++)
+          {
+            std::complex<T>& c = new_data[alms.index(l,m)];
+            c.real() = gsl_ran_gaussian(rng, Al); 
+            c.imag() = gsl_ran_gaussian(rng, Al); 
+          }
+      }
+    return r;
+  } 
+};
+
+#endif

src/loadGadget.cpp

@@ -256,7 +256,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i
       }
     catch (const InvalidUnformattedAccess& e)
       {
-	cerr << "Invalid formatted while reading ID" << endl;
+	cerr << "Invalid unformatted access while reading ID" << endl;
 	delete f;
 	delete data;
 	return 0;

src/yorick_nc3.cpp

@@ -1,5 +1,11 @@
+#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>

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);
+  
+}