From 86dc4bd2493060e26c3cc169a7383521f0dadf22 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Thu, 27 Feb 2014 15:30:21 +0100
Subject: [PATCH 01/19] Fixed default arguments in healpix_transform. Fixed
 fftw detection in CMake

---
 CMakeLists.txt                            | 4 ++++
 sample/CMakeLists.txt                     | 8 +++++++-
 src/fourier/details/healpix_transform.hpp | 2 +-
 3 files changed, 12 insertions(+), 2 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index c4888a4..7160419 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -61,6 +61,10 @@ pkg_check_modules(FFTW3 fftw3>=3.3)
 pkg_check_modules(FFTW3F fftw3f>=3.3)
 pkg_check_modules(EIGEN3 eigen3)
 
+
+find_library(FFTW3F_LIBRARY_FULL fftw3f PATHS ${FFTW3F_LIBDIR} NO_DEFAULT_PATH)
+find_library(FFTW3_LIBRARY_FULL fftw3 PATHS ${FFTW3_LIBDIR} NO_DEFAULT_PATH)
+
 include(FindPackageHandleStandardArgs)
 set(NETCDF_FIND_REQUIRED TRUE)
 set(GSL_FIND_REQUIRED TRUE)
diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt
index 969480c..64bb51f 100644
--- a/sample/CMakeLists.txt
+++ b/sample/CMakeLists.txt
@@ -48,8 +48,14 @@ add_executable(testBSP testBSP.cpp)
 target_link_libraries(testBSP ${tolink})
 
 if (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND)
+  IF (FFTW3F_LIBRARY_FULL)
+    SET(FFTW3_LIB ${FFTW3F_LIBRARY_FULL})
+  ENDIF (FFTW3F_LIBRARY_FULL)
+  IF (FFTW3_LIBRARY_FULL)
+    SET(FFTW3_LIB ${FFTW3_LIB} ${FFTW3_LIBRARY_FULL})
+  ENDIF (FFTW3_LIBRARY_FULL)
   add_executable(test_fft_calls test_fft_calls.cpp)
-  target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES})
+  target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIB})
 endif (FFTW3_FOUND AND FFTW3F_FOUND AND  EIGEN3_FOUND)
 
 if (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
diff --git a/src/fourier/details/healpix_transform.hpp b/src/fourier/details/healpix_transform.hpp
index d3f7927..37e79e8 100644
--- a/src/fourier/details/healpix_transform.hpp
+++ b/src/fourier/details/healpix_transform.hpp
@@ -67,7 +67,7 @@ namespace CosmoTool
       sharp_make_triangular_alm_info (Lmax, Mmax, 1, &ainfo);
     }
 
-    HealpixFourierTransform(long nSide, long Lmax, long Mmax, int iterate = 0, const std::valarray<double>& weights )
+    HealpixFourierTransform(long nSide, long Lmax, long Mmax, int iterate, const std::valarray<double>& weights )
       : realMap(nSide), fourierMap(Lmax, Mmax), ainfo(0), ginfo(0), m_iterate(iterate)
     {
       sharp_make_weighted_healpix_geom_info (nSide, 1, &weights[0], &ginfo);

From aed2920be7b28532b567bebe2ae376f136b2b7fa Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Thu, 27 Feb 2014 15:59:25 +0100
Subject: [PATCH 02/19] Fixed interface get_Kvec_p

---
 sample/test_fft_calls.cpp              | 1 +
 src/fourier/details/euclidian_maps.hpp | 4 ++--
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/sample/test_fft_calls.cpp b/sample/test_fft_calls.cpp
index b98fcea..e49ad82 100644
--- a/sample/test_fft_calls.cpp
+++ b/sample/test_fft_calls.cpp
@@ -89,5 +89,6 @@ int main(int argc, char **argv)
   test_2d<float>(128,128);
   test_2d<float>(128,131);  
   test_2d<float>(128,130);  
+  test_2d<float>(131,128);  
   return 0;
 }
diff --git a/src/fourier/details/euclidian_maps.hpp b/src/fourier/details/euclidian_maps.hpp
index 53b0278..5b9b95d 100644
--- a/src/fourier/details/euclidian_maps.hpp
+++ b/src/fourier/details/euclidian_maps.hpp
@@ -172,7 +172,7 @@ namespace CosmoTool
     }
 
     template<typename Array2>
-    void get_Kvec(long p, Array2& kvec) const
+    void get_Kvec_p(long p, Array2& kvec) const
     {
       const DimArray& dims = this->getDims();
       DimArray d(delta_k.size());
@@ -212,7 +212,7 @@ namespace CosmoTool
       return std::sqrt(k2);
     }
 
-    double get_K(long p) const
+    double get_K_p(long p) const
     {
       const DimArray& dims = this->getDims();
       DimArray d(delta_k.size());

From c72470024668c4defdc8701231e8f388f79887ef Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Thu, 27 Feb 2014 16:03:42 +0100
Subject: [PATCH 03/19] Fixed calls to get_K_p

---
 src/fourier/details/euclidian_spectrum_1d.hpp | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/fourier/details/euclidian_spectrum_1d.hpp b/src/fourier/details/euclidian_spectrum_1d.hpp
index 36b524a..4d6fb9b 100644
--- a/src/fourier/details/euclidian_spectrum_1d.hpp
+++ b/src/fourier/details/euclidian_spectrum_1d.hpp
@@ -115,7 +115,7 @@ namespace CosmoTool
 
     for (long p = 1; p < rand_map.size(); p++)
       {
-	double A_k = std::sqrt(0.5*V*f(rand_map.get_K(p)));
+	double A_k = std::sqrt(0.5*V*f(rand_map.get_K_p(p)));
 	d[p] = std::complex<T>(gsl_ran_gaussian(rng, A_k),
 			       gsl_ran_gaussian(rng, A_k));
       }
@@ -164,7 +164,7 @@ namespace CosmoTool
     std::complex<T> *d = m.data();
 
     for (long p = 0; p < m_c.size(); p++)
-      d[p] *= f(m_c.get_K(p));
+      d[p] *= f(m_c.get_K_p(p));
   }
 
   template<typename T>
@@ -174,7 +174,7 @@ namespace CosmoTool
     std::complex<T> *d = m.data();
 
     for (long p = 0; p < m_c.size(); p++)
-      d[p] *= std::sqrt(f(m_c.get_K(p)));
+      d[p] *= std::sqrt(f(m_c.get_K_p(p)));
   }
 
   template<typename T>
@@ -185,7 +185,7 @@ namespace CosmoTool
 
     for (long p = 0; p < m_c.size(); p++)
      {
-        T A = f(m_c.get_K(p));
+        T A = f(m_c.get_K_p(p));
         if (A==0)
           d[p] = 0;
         else
@@ -201,7 +201,7 @@ namespace CosmoTool
 
     for (long p = 0; p < m_c.size(); p++)
       {
-        T A = std::sqrt(f(m_c.get_K(p)));
+        T A = std::sqrt(f(m_c.get_K_p(p)));
         if (A == 0)
           d[p] = 0;
         else

From 126e1461b9f6362ae7e820dbd0a9ad208477e22c Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Mon, 3 Mar 2014 18:27:38 +0100
Subject: [PATCH 04/19] Attempts to fix the conjugation in the dot-product

---
 src/cosmopower.cpp                            | 2 +-
 src/fourier/details/euclidian_maps.hpp        | 4 ++--
 src/fourier/details/euclidian_spectrum_1d.hpp | 9 ++++++++-
 3 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp
index a96a6bb..20026a2 100644
--- a/src/cosmopower.cpp
+++ b/src/cosmopower.cpp
@@ -212,7 +212,7 @@ double CosmoPower::powerBDM(double k)
 
 double CosmoPower::powerTest(double k)
 {
-  return 1/(1+k*k);
+  return normPower;//1/(1+k*k);
 }
 
 /*
diff --git a/src/fourier/details/euclidian_maps.hpp b/src/fourier/details/euclidian_maps.hpp
index 5b9b95d..5673565 100644
--- a/src/fourier/details/euclidian_maps.hpp
+++ b/src/fourier/details/euclidian_maps.hpp
@@ -254,8 +254,8 @@ namespace CosmoTool
 	  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];
+	      result += 2*conj(d1[q0]) * d2[q0];
+	      result += 2*conj(d1[q1]) * d2[q1];
 	    }
 	}
       return result;
diff --git a/src/fourier/details/euclidian_spectrum_1d.hpp b/src/fourier/details/euclidian_spectrum_1d.hpp
index 4d6fb9b..01a7dac 100644
--- a/src/fourier/details/euclidian_spectrum_1d.hpp
+++ b/src/fourier/details/euclidian_spectrum_1d.hpp
@@ -138,7 +138,7 @@ namespace CosmoTool
 	plane_size *= dims[q];
       }
 
-    for (long p = 1; p < plane_size/2; p++)
+    for (long p = 1; p < plane_size/2+1; p++)
       {
 	long q = (p+1)*dims[0]-1;
 	long q2 = (plane_size-p+1)*dims[0]-1;
@@ -147,6 +147,13 @@ namespace CosmoTool
 	d[q] = conj(d[q2]);
       }
 
+    for (long p = 1; p < plane_size/2+1; p++)
+      {
+	long q = (p)*dims[0];
+	long q2 = (plane_size-p)*dims[0];
+	d[q] = conj(d[q2]);
+      }
+
     if (alleven)
       {
 	long q = 0;

From a2373d8919f5dea4c571824cf3e59ce9d3b9d1d6 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Mon, 3 Mar 2014 19:06:31 +0100
Subject: [PATCH 05/19] Fixed types

---
 src/fourier/details/euclidian_maps.hpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/fourier/details/euclidian_maps.hpp b/src/fourier/details/euclidian_maps.hpp
index 5673565..4cbf6a3 100644
--- a/src/fourier/details/euclidian_maps.hpp
+++ b/src/fourier/details/euclidian_maps.hpp
@@ -36,6 +36,7 @@ knowledge of the CeCILL license and that you accept its terms.
 #ifndef __DETAILS_EUCLIDIAN_MAPS
 #define __DETAILS_EUCLIDIAN_MAPS
 
+#include <cmath>
 
 namespace CosmoTool
 {
@@ -246,7 +247,7 @@ namespace CosmoTool
 	    {
 	      long idx = q0+dims[0]*p;
 	      assert(idx < this->size());
-	      result += 2*(conj(d1[idx]) * d2[idx]).real();
+	      result += T(2)*(std::conj(d1[idx]) * d2[idx]).real();
 	    }
 	}
       if (even0)
@@ -254,8 +255,8 @@ namespace CosmoTool
 	  for (long p = 0; p < plane_size; p++)
 	    {
 	      long q0 = N0*p, q1 = (p+1)*N0-1;
-	      result += 2*conj(d1[q0]) * d2[q0];
-	      result += 2*conj(d1[q1]) * d2[q1];
+	      result += T(2)*std::conj(d1[q0]) * d2[q0];
+	      result += T(2)*std::conj(d1[q1]) * d2[q1];
 	    }
 	}
       return result;

From 9ae42a2442f0dccf89422e77baf0725749645a43 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Wed, 2 Apr 2014 09:52:06 +0200
Subject: [PATCH 06/19] Fixed unit_vel in Ramses. Added a new type alias in
 FourierMap

---
 src/fourier/base_types.hpp | 3 +++
 src/loadRamses.cpp         | 2 +-
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/fourier/base_types.hpp b/src/fourier/base_types.hpp
index ff092ef..42b579a 100644
--- a/src/fourier/base_types.hpp
+++ b/src/fourier/base_types.hpp
@@ -55,6 +55,7 @@ namespace CosmoTool
   protected:
     SpectrumFunction() {}
   public:
+    typedef T type;
     typedef Eigen::Array<T, 1, Eigen::Dynamic> VecType;
     typedef Eigen::Map<VecType, Eigen::Aligned> MapType;
     typedef Eigen::Map<VecType const, Eigen::Aligned> ConstMapType;
@@ -88,6 +89,7 @@ namespace CosmoTool
     FourierMap() {}
     
   public:
+    typedef T type;
     typedef Eigen::Array<T, 1, Eigen::Dynamic> VecType;
     typedef Eigen::Map<VecType, Eigen::Aligned> MapType;
     typedef Eigen::Map<VecType const, Eigen::Aligned> ConstMapType;
@@ -190,6 +192,7 @@ namespace CosmoTool
   class MapUtilityFunction
   {
   public:
+    typedef T type;
     typedef SpectrumFunction<T> Spectrum;
     typedef boost::shared_ptr<Spectrum> Spectrum_ptr;
     typedef FourierMap<std::complex<T> > FMap;
diff --git a/src/loadRamses.cpp b/src/loadRamses.cpp
index 3eb909c..5154921 100644
--- a/src/loadRamses.cpp
+++ b/src/loadRamses.cpp
@@ -309,7 +309,7 @@ CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputI
 
   double hubble = info.aexp*info.aexp/info.unit_t / (1e5/CM_IN_MPC);
   double L0 = info.boxSize*info.unitLength*hubble/100/CM_IN_MPC/info.aexp;
-  double unit_vel = L0*hubble/info.aexp;
+  double unit_vel = L0*100/info.aexp;
   
   while (1)
     {

From bba009ef42489c6d4fa71b5aa19108545716ca16 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Wed, 2 Apr 2014 16:53:14 +0200
Subject: [PATCH 07/19] Export correctly the const-ness of saveArray

---
 src/yorick.hpp     | 2 +-
 src/yorick_nc3.cpp | 8 ++++----
 src/yorick_nc4.cpp | 8 ++++----
 3 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/yorick.hpp b/src/yorick.hpp
index 39fd20d..59d4656 100644
--- a/src/yorick.hpp
+++ b/src/yorick.hpp
@@ -218,7 +218,7 @@ namespace CosmoTool
 
   template<typename T>
   void saveArray(const std::string& fname,
-		 T *array, uint32_t *dimList, uint32_t rank);
+		 const T *array, uint32_t *dimList, uint32_t rank);
 
   template<typename T>
   void loadArray(const std::string& fname,
diff --git a/src/yorick_nc3.cpp b/src/yorick_nc3.cpp
index d80781e..8a69d9b 100644
--- a/src/yorick_nc3.cpp
+++ b/src/yorick_nc3.cpp
@@ -235,7 +235,7 @@ namespace CosmoTool {
 
   template<typename T>
   void saveArray(const std::string& fname,
-		 T *array, uint32_t *dimList, uint32_t rank)
+		 const T *array, uint32_t *dimList, uint32_t rank)
   {
     NcFile f(fname.c_str(), NcFile::Replace, 0, 0, NcFile::Netcdf4);
     
@@ -300,10 +300,10 @@ namespace CosmoTool {
 				  double*& array, uint32_t *&dimList, uint32_t& rank);
 
   template void saveArray<int>(const std::string& fname,
-			       int *array, uint32_t *dimList, uint32_t rank);
+			       const int *array, uint32_t *dimList, uint32_t rank);
   template void saveArray<float>(const std::string& fname,
-				 float *array, uint32_t *dimList, uint32_t rank);
+				 const float *array, uint32_t *dimList, uint32_t rank);
   template void saveArray<double>(const std::string& fname,
-				  double *array, uint32_t *dimList, uint32_t rank);
+				  const double *array, uint32_t *dimList, uint32_t rank);
   
 }
diff --git a/src/yorick_nc4.cpp b/src/yorick_nc4.cpp
index ed4cd9e..49ae5b5 100644
--- a/src/yorick_nc4.cpp
+++ b/src/yorick_nc4.cpp
@@ -203,7 +203,7 @@ namespace CosmoTool {
 
   template<typename T>
   void saveArray(const std::string& fname,
-		 T *array, uint32_t *dimList, uint32_t rank)
+		 const T *array, uint32_t *dimList, uint32_t rank)
   {
     NcFile f(fname.c_str(), NcFile::replace);
     
@@ -263,10 +263,10 @@ namespace CosmoTool {
 				  double*& array, uint32_t *&dimList, uint32_t& rank);
 
   template void saveArray<int>(const std::string& fname,
-			       int *array, uint32_t *dimList, uint32_t rank);
+			       const int *array, uint32_t *dimList, uint32_t rank);
   template void saveArray<float>(const std::string& fname,
-				 float *array, uint32_t *dimList, uint32_t rank);
+				 const float *array, uint32_t *dimList, uint32_t rank);
   template void saveArray<double>(const std::string& fname,
-				  double *array, uint32_t *dimList, uint32_t rank);
+				  const double *array, uint32_t *dimList, uint32_t rank);
   
 }

From e3ff41e30ae1de71bb42d09c60f5a5228de6f4d8 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Wed, 2 Apr 2014 17:56:09 +0200
Subject: [PATCH 08/19] Added new calls specific to EuclidianTransform

---
 src/fourier/details/euclidian_transform.hpp | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/src/fourier/details/euclidian_transform.hpp b/src/fourier/details/euclidian_transform.hpp
index eaab4fd..e5ef34a 100644
--- a/src/fourier/details/euclidian_transform.hpp
+++ b/src/fourier/details/euclidian_transform.hpp
@@ -122,12 +122,22 @@ namespace CosmoTool
       calls::execute(m_synthesis);
       realMap->scale(1/volume);
     }
+
+    void synthesis_unnormed()
+    {
+      calls::execute(m_synthesis);
+    }
     
     void analysis()
     {
       calls::execute(m_analysis);
       fourierMap->scale(volume/N);
     }
+
+    void analysis_unnormed()
+    {
+      calls::execute(m_analysis);
+    }
     
     void synthesis_conjugate()
     {

From b80734bf07ad698b1daed93f42c461d20418c928 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Thu, 10 Apr 2014 09:12:01 +0200
Subject: [PATCH 09/19] Enforce C compiler as given by CMake

---
 CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 858dbaa..9d69ebe 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -37,7 +37,7 @@ if (ENABLE_SHARP)
   ExternalProject_Add(sharp
     SOURCE_DIR ${SHARP_SOURCE}
     BUILD_IN_SOURCE 1 
-    CONFIGURE_COMMAND ${SHARP_SOURCE}/configure --prefix=${DEP_BUILD}
+    CONFIGURE_COMMAND ${SHARP_SOURCE}/configure "CC=${CMAKE_C_COMPILER}" "CXX=${CMAKE_CXX_COMPILER}" --prefix=${DEP_BUILD}
     BUILD_COMMAND ${CMAKE_MAKE_PROGRAM}
     INSTALL_COMMAND echo "No install"
   )

From cc769b2b1d525347edcfebe3394a76b0a6631a68 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Mon, 5 May 2014 12:07:55 +0200
Subject: [PATCH 10/19] Added state definition in SPHSmooth. Fixed cmake
 LINK_FLAGS problem

---
 CMakeLists.txt        |  12 +++++
 sample/CMakeLists.txt |   2 +-
 sample/testSmooth.cpp |   6 ++-
 src/sphSmooth.hpp     |  43 ++++++++-------
 src/sphSmooth.tcc     | 120 +++++++++++++++++++++++-------------------
 5 files changed, 106 insertions(+), 77 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9d69ebe..e2e45f6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -71,6 +71,18 @@ set(GSL_FIND_REQUIRED TRUE)
 FIND_PACKAGE_HANDLE_STANDARD_ARGS(NetCDF DEFAULT_MSG NETCDF_LIBRARY NETCDFCPP_LIBRARY NETCDF_INCLUDE_PATH)
 FIND_PACKAGE_HANDLE_STANDARD_ARGS(GSL DEFAULT_MSG GSL_LIBRARY GSLCBLAS_LIBRARY GSL_INCLUDE_PATH)
 
+IF(ENABLE_OPENMP)
+
+  IF (NOT OPENMP_FOUND)
+    MESSAGE(ERROR "No known compiler option for enabling OpenMP")
+  ENDIF(NOT OPENMP_FOUND)
+
+  SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+  SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+  SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKED_FLAGS} ${OpenMP_C_FLAGS}")
+ENDIF(ENABLE_OPENMP)
+
+
 
 set(GSL_LIBRARIES ${GSL_LIBRARY} ${GSLCBLAS_LIBRARY})
 
diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt
index 64bb51f..5b58290 100644
--- a/sample/CMakeLists.txt
+++ b/sample/CMakeLists.txt
@@ -62,7 +62,7 @@ if (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
   include_directories(${SHARP_INCLUDE_PATH})
   add_executable(test_healpix_calls test_healpix_calls.cpp)
   target_link_libraries(test_healpix_calls ${tolink} ${SHARP_LIBRARIES})
-  set_target_properties(test_healpix_calls PROPERTIES COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS})
+  set_target_properties(test_healpix_calls PROPERTIES COMPILE_FLAGS "${OpenMP_CXX_FLAGS}" LINK_FLAGS "${OpenMP_CXX_FLAGS}")
   add_dependencies(test_healpix_calls sharp)
 endif (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
 
diff --git a/sample/testSmooth.cpp b/sample/testSmooth.cpp
index 4f3011d..ddb8c3e 100644
--- a/sample/testSmooth.cpp
+++ b/sample/testSmooth.cpp
@@ -85,14 +85,16 @@ int main()
   uint32_t dims[] = { NX, NX };
   ProgressiveOutput<ComputePrecision> out = 
     ProgressiveOutput<ComputePrecision>::saveArrayProgressive("out.nc", dims, 2);
+//#pragma omp parallel for schedule(static)
   for (uint32_t iy = 0; iy < NX; iy++)
     {
+      MySmooth::SPHState state;
       cout << "iy=" << iy << endl;
       for (uint32_t ix = 0; ix < NX; ix++)
         {
            MyTree::coords c = { 1.0*ix/NX, 1.0*iy/NX };
-           smooth.fetchNeighbours(c);
-           out.put(smooth.computeSmoothedValue(c, unit_fun));
+           smooth.fetchNeighbours(c, &state);
+           out.put(smooth.computeSmoothedValue(c, unit_fun, &state));
 
 	  
         }
diff --git a/src/sphSmooth.hpp b/src/sphSmooth.hpp
index f6fe582..b347178 100644
--- a/src/sphSmooth.hpp
+++ b/src/sphSmooth.hpp
@@ -35,7 +35,7 @@ knowledge of the CeCILL license and that you accept its terms.
 
 #ifndef __COSMOTOOL_SPH_SMOOTH_HPP
 #define __COSMOTOOL_SPH_SMOOTH_HPP
-
+#include <boost/shared_ptr.hpp>
 #include "config.hpp"
 #include "mykdtree.hpp"
 
@@ -60,37 +60,48 @@ namespace CosmoTool
     typedef void (*runParticleValue)(ValType& t);
 
   public:
+    struct SPHState
+    {
+      boost::shared_ptr<SPHCell *[]> ngb;
+      boost::shared_ptr<CoordType[]> distances;
+      typename SPHTree::coords currentCenter;
+      int currentNgb;
+      ComputePrecision smoothRadius;
+    };
+    
+  
     SPHSmooth(SPHTree *tree, uint32_t Nsph);
     virtual ~SPHSmooth();   
 
-    void fetchNeighbours(const typename SPHTree::coords& c);
+    void fetchNeighbours(const typename SPHTree::coords& c, SPHState *state = 0);
+    
     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;
+      return internal.currentCenter;
     }
 
     ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c,
-					  computeParticleValue fun);    
+					  computeParticleValue fun, SPHState *state = 0);    
     ComputePrecision computeInterpolatedValue(const typename SPHTree::coords& c,
-					      computeParticleValue fun);    
+					      computeParticleValue fun, SPHState *state = 0);    
     ComputePrecision getMaxDistance(const typename SPHTree::coords& c, 
 				    SPHNode *node) const;
 
     ComputePrecision getSmoothingLen() const
     {
-      return smoothRadius;
+      return internal.smoothRadius;
     }
 
     // TO USE WITH EXTREME CARE !
     void setSmoothingLen(ComputePrecision len)
     {
-      smoothRadius = len;
+      internal.smoothRadius = len;
     }
     // END
 
-    void runForEachNeighbour(runParticleValue fun);
+    void runForEachNeighbour(runParticleValue fun, SPHState *state = 0);
     void addGridSite(const typename SPHTree::coords& c);
 
     bool hasNeighbours() const;
@@ -100,32 +111,26 @@ namespace CosmoTool
     SPHTree *getTree() { return tree; }
 
     void changeNgb(uint32_t newMax) {
-      delete[] ngb;
-      delete[] distances;
-      ngb = new SPHCell *[newMax];
-      distances = new CoordType[newMax];
+      internal.ngb = boost::shared_ptr<SPHCell *>(new SPHCell *[newMax]);
+      internal.distances = boost::shared_ptr<CoordType>(new CoordType[newMax]);
       maxNgb = newMax;
     }
 
-    uint32_t getCurrent() const { return currentNgb; }
+    uint32_t getCurrent() const { return internal.currentNgb; }
 
     uint32_t getNgb() const { return maxNgb; }
  
   protected:
-    SPHCell **ngb;
-    CoordType *distances;
+    SPHState internal;
     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);
+				   computeParticleValue fun, SPHState *state);
     void runUnrollNode(SPHNode *node,
 		       runParticleValue fun);
   };
diff --git a/src/sphSmooth.tcc b/src/sphSmooth.tcc
index 0cf7de8..b9c9d2e 100644
--- a/src/sphSmooth.tcc
+++ b/src/sphSmooth.tcc
@@ -1,4 +1,5 @@
 #include <cmath>
+#include "algo.hpp"
 
 namespace CosmoTool {
 
@@ -7,29 +8,27 @@ SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph)
 {
   this->Nsph = Nsph;
   this->tree = tree;
-  this->currentNgb = 0;
+  internal.currentNgb = 0;
 
   this->maxNgb = Nsph;
-  ngb = new SPHCell *[maxNgb];  
-  distances = new CoordType[maxNgb];
+  internal.ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[maxNgb]);  
+  internal.distances = boost::shared_ptr<CoordType[]>(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)
+							 computeParticleValue fun, SPHState *state)
 {
   CoordType weight;
 
-  d /= smoothRadius;
+  d /= state->smoothRadius;
   weight = getKernel(d);
 
   if (cell.val.weight != 0)
@@ -47,86 +46,91 @@ SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, uin
 
   if (requested > maxNgb)
     {
-      delete[] ngb;
-      delete[] distances;
-
       maxNgb = requested;
-      ngb = new SPHCell *[maxNgb];
-      distances = new CoordType[maxNgb];
+      internal.ngb = new SPHCell *[maxNgb];
+      internal.distances = new CoordType[maxNgb];
     }
 
-  memcpy(currentCenter, c, sizeof(c));
-  tree->getNearestNeighbours(c, requested, ngb, distances);
+  memcpy(internal.currentCenter, c, sizeof(c));
+  tree->getNearestNeighbours(c, requested, internal.ngb, internal.distances);
 
-  currentNgb = 0;
-  for (uint32_t i = 0; i < requested && ngb[i] != 0; i++,currentNgb++)
+  internal.currentNgb = 0;
+  for (uint32_t i = 0; i < requested && internal.ngb[i] != 0; i++,internal.currentNgb++)
     {
-      distances[i] = sqrt(distances[i]);
-      d2 = distances[i];
+      internal.distances[i] = sqrt(internal.distances[i]);
+      d2 = internal.distances[i];
       if (d2 > max_dist)
         max_dist = d2;
     }
 
-  smoothRadius = max_dist / 2;  
+  internal.smoothRadius = max_dist / 2;  
 }
 
 template<typename ValType, int Ndims>
-void
-SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c)
+void SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, SPHState *state)
 {
   ComputePrecision d2, max_dist = 0;
   uint32_t requested = Nsph;
 
-  memcpy(currentCenter, c, sizeof(c));
-  tree->getNearestNeighbours(c, requested, ngb, distances);
+  if (state != 0) {
+    state->distances = boost::shared_ptr<CoordType[]>(new CoordType[Nsph]);
+    state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[Nsph]);
+  } else
+    state = &internal;
+    
+  memcpy(state->currentCenter, c, sizeof(c));
+  
+  tree->getNearestNeighbours(c, requested, state->ngb.get(), state->distances.get());
 
-  currentNgb = 0;
-  for (uint32_t i = 0; i < requested && ngb[i] != 0; i++,currentNgb++)
+  state->currentNgb = 0;
+  for (uint32_t i = 0; i < requested && state->ngb[i] != 0; i++,state->currentNgb++)
     {
-      distances[i] = sqrt(distances[i]);
-      d2 = distances[i];
+      d2 = internal.distances[i] = sqrt(internal.distances[i]);
       if (d2 > max_dist)
-	max_dist = d2;
+        max_dist = d2;
     }
 
-  smoothRadius = max_dist / 2;  
-}					
+  state->smoothRadius = max_dist / 2;  
+}
+
 
 template<typename ValType, int Ndims>
 void
 SPHSmooth<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords& c,
-						  ComputePrecision radius)
+                                                  ComputePrecision radius)
 {
   uint32_t numPart;
   ComputePrecision d2, max_dist = 0;
 
-  memcpy(currentCenter, c, sizeof(c));
+  memcpy(internal.currentCenter, c, sizeof(c));
 
-  currentNgb = tree->getIntersection(c, radius, ngb, distances,
+  internal.currentNgb = tree->getIntersection(c, radius, internal.ngb, internal.distances,
 				  maxNgb);
 
-  for (uint32_t i = 0; i < currentNgb; i++)
+  for (uint32_t i = 0; i < internal.currentNgb; i++)
     {
-      distances[i] = sqrt(distances[i]);
-      d2 = distances[i];
+      d2 = internal.distances[i] = sqrt(internal.distances[i]);
       if (d2 > max_dist)
-	max_dist = d2;
+        max_dist = d2;
     }
-  smoothRadius = max_dist / 2;  
+  internal.smoothRadius = max_dist / 2;  
 }				
 
 template<typename ValType, int Ndims>
 ComputePrecision 
 SPHSmooth<ValType,Ndims>::computeSmoothedValue(const typename SPHTree::coords& c,
-					       computeParticleValue fun)
+                                               computeParticleValue fun, SPHState *state)
 {
+  if (state == 0)
+    state = &internal;
+    
   ComputePrecision outputValue = 0;
   ComputePrecision max_dist = 0;
-  ComputePrecision r3 = smoothRadius * smoothRadius * smoothRadius;
+  ComputePrecision r3 = cube(state->smoothRadius);
 
-  for (uint32_t i = 0; i < currentNgb; i++)
+  for (uint32_t i = 0; i < state->currentNgb; i++)
     {
-      outputValue += computeWValue(c, *ngb[i], distances[i], fun);
+      outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun, state);
     }
 
   return outputValue / r3;
@@ -141,27 +145,33 @@ ComputePrecision interpolateOne(const ValType& t)
 // WARNING ! Cell's weight must be 1 !!!
 template<typename ValType, int Ndims>
 ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typename SPHTree::coords& c,
-								    computeParticleValue fun)
+								    computeParticleValue fun, SPHState *state)
 {
+  if (state == 0)
+    state = &internal;
+    
   ComputePrecision outputValue = 0;
   ComputePrecision max_dist = 0;
   ComputePrecision weight = 0;
 
-  for (uint32_t i = 0; i < currentNgb; i++)
+  for (uint32_t i = 0; i < state->currentNgb; i++)
     {
-      outputValue += computeWValue(c, *ngb[i], distances[i], fun);
-      weight += computeWValue(c, *ngb[i], distances[i], interpolateOne);
+      outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun);
+      weight += computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne);
     }
 
   return (outputValue == 0) ? 0 : (outputValue / weight);  
 }
 
 template<typename ValType, int Ndims>
-void SPHSmooth<ValType,Ndims>::runForEachNeighbour(runParticleValue fun)
+void SPHSmooth<ValType,Ndims>::runForEachNeighbour(runParticleValue fun, SPHState *state)
 {
-  for (uint32_t i = 0; i < currentNgb; i++)
+  if (state == 0)
+    state = &internal;
+    
+  for (uint32_t i = 0; i < state->currentNgb; i++)
     {
-      fun(ngb[i]);
+      fun(state->ngb[i]);
     }
 }
 
@@ -172,13 +182,13 @@ void SPHSmooth<ValType,Ndims>::addGridSite(const typename SPHTree::coords& c)
   ComputePrecision outputValue = 0;
   ComputePrecision max_dist = 0;
   
-  ComputePrecision r3 = smoothRadius * smoothRadius * smoothRadius;
+  ComputePrecision r3 = cube(internal.smoothRadius);
 
-  for (uint32_t i = 0; i < currentNgb; i++)
+  for (uint32_t i = 0; i < internal.currentNgb; i++)
     {
-      ComputePrecision d = distances[i];
-      SPHCell& cell = *ngb[i];
-      cell.val.weight += getKernel(d/smoothRadius) / r3;
+      ComputePrecision d = internal.distances[i];
+      SPHCell& cell = *internal.ngb[i];
+      cell.val.weight += getKernel(d/internal.smoothRadius) / r3;
 
     }
 }
@@ -202,7 +212,7 @@ SPHSmooth<ValType,Ndims>::getKernel(ComputePrecision x) const
 template<typename ValType, int Ndims>
 bool SPHSmooth<ValType,Ndims>::hasNeighbours() const
 {
-  return (currentNgb != 0);
+  return (internal.currentNgb != 0);
 }
 
 template<class ValType1, class ValType2, int Ndims>

From 774dc99e67a14c097bfca373c57162eab845181b Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Mon, 5 May 2014 12:12:09 +0200
Subject: [PATCH 11/19] Temporarily disable testSmooth as it depends on Boost

---
 sample/CMakeLists.txt | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt
index 5b58290..2053060 100644
--- a/sample/CMakeLists.txt
+++ b/sample/CMakeLists.txt
@@ -11,8 +11,8 @@ target_link_libraries(testBQueue ${tolink})
 add_executable(testInterpolate testInterpolate.cpp)
 target_link_libraries(testInterpolate ${tolink})
 
-add_executable(testSmooth testSmooth.cpp)
-target_link_libraries(testSmooth ${tolink})
+#add_executable(testSmooth testSmooth.cpp)
+#target_link_libraries(testSmooth ${tolink})
 
 add_executable(testkd testkd.cpp)
 target_link_libraries(testkd ${tolink})

From 4a7859587ee757c018c63ec6b039b62fbe1e46a9 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem.lavaux@laposte.net>
Date: Tue, 20 May 2014 09:25:32 +0200
Subject: [PATCH 12/19] Imported the simple3DFilter to the sample directory

---
 sample/CMakeLists.txt     |   5 +-
 sample/simple3DFilter.cpp | 193 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 197 insertions(+), 1 deletion(-)
 create mode 100644 sample/simple3DFilter.cpp

diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt
index 64bb51f..9475f13 100644
--- a/sample/CMakeLists.txt
+++ b/sample/CMakeLists.txt
@@ -34,7 +34,7 @@ target_link_libraries(testPool ${tolink})
 
 if (HDF5_FOUND)
 	  add_executable(testReadFlash testReadFlash.cpp)
-	  target_link_libraries(testReadFlash ${tolink})
+    target_link_libraries(testReadFlash ${tolink})
 endif (HDF5_FOUND)
 
 
@@ -69,3 +69,6 @@ endif (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
 add_executable(test_cosmopower test_cosmopower.cpp)
 target_link_libraries(test_cosmopower ${tolink})
 
+
+add_executable(simple3DFilter simple3DFilter.cpp)
+target_link_libraries(simple3DFilter ${tolink})
diff --git a/sample/simple3DFilter.cpp b/sample/simple3DFilter.cpp
new file mode 100644
index 0000000..8b91332
--- /dev/null
+++ b/sample/simple3DFilter.cpp
@@ -0,0 +1,193 @@
+#include <cassert>
+#include <CosmoTool/yorick.hpp>
+#include <CosmoTool/sphSmooth.hpp>
+#include <CosmoTool/mykdtree.hpp>
+#include <CosmoTool/miniargs.hpp>
+
+using namespace std;
+using namespace CosmoTool;
+
+#define N_SPH 16
+
+struct VCoord{
+  float v[3];
+};
+
+template<int i>
+ComputePrecision getVelocity(const VCoord& v)
+{
+  return v.v[i];
+}
+
+ComputePrecision getUnity(const VCoord& v)
+{
+  return 1.0;
+}
+
+typedef SPHSmooth<VCoord> MySmooth;
+typedef MySmooth::SPHTree MyTree;
+typedef MyTree::Cell MyCell;
+
+int main(int argc, char **argv)
+{
+  char *fname1, *fname2;
+  double rLimit, boxsize, rLimit2, cx, cy, cz;
+  int Nres;
+
+  MiniArgDesc args[] = {
+    { "INPUT DATA1", &fname1, MINIARG_STRING },
+    { "RADIUS LIMIT", &rLimit, MINIARG_DOUBLE },
+    { "BOXSIZE", &boxsize, MINIARG_DOUBLE },
+    { "RESOLUTION", &Nres, MINIARG_INT },
+    { "CX", &cx, MINIARG_DOUBLE },
+    { "CY", &cy, MINIARG_DOUBLE },
+    { "CZ", &cz, MINIARG_DOUBLE },
+    { 0, 0, MINIARG_NULL }
+  };
+
+  if (!parseMiniArgs(argc, argv, args))
+    return 1;
+
+  float *v1_data;
+  uint32_t *dimList;
+  uint32_t rank;
+  uint32_t N1_points, N2_points;
+  int *bins = new int[Nres*Nres*Nres];
+
+  rLimit2 = rLimit*rLimit;
+
+  loadArray(fname1, v1_data, dimList, rank);
+  assert(rank == 2);
+  assert(dimList[1] == 6);
+
+  N1_points = dimList[0];
+  delete[] dimList;
+  
+  cout << "Got " << N1_points << " in the first file." << endl;
+
+  MyCell *allCells_1 = new MyCell[N1_points];
+  
+  for (long i = 0; i < Nres*Nres*Nres; i++)
+    bins[i] = 0;
+
+  cout << "Shuffling data in cells..." << endl;
+  for (int i = 0 ; i < N1_points; i++)
+    {
+      for (int j = 0; j < 3; j++)
+	allCells_1[i].coord[j] = v1_data[i*6 + j];
+      for (int k = 0; k < 3; k++)
+	allCells_1[i].val.pValue.v[k] = v1_data[i*6 + 3 + k];
+      allCells_1[i].active = true;
+      allCells_1[i].val.weight = 0.0;
+
+      long rx = floor((allCells_1[i].coord[0]+cx)*Nres/boxsize+0.5);
+      long ry = floor((allCells_1[i].coord[1]+cy)*Nres/boxsize+0.5);
+      long rz = floor((allCells_1[i].coord[2]+cz)*Nres/boxsize+0.5);
+      
+      if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres)
+	  continue;
+	  
+      bins[rx + ry*Nres + rz*Nres*Nres]++;
+    }
+  delete[] v1_data;
+   uint32_t dims[3] = { Nres, Nres, Nres } ;
+   saveArray("num_in_cell.nc", bins, dims, 3);
+
+  cout << "Building trees..." << endl;
+  MyTree tree1(allCells_1, N1_points);
+
+  cout << "Creating smoothing filter..." << endl;
+  MySmooth smooth1(&tree1, N_SPH);
+
+  uint32_t outDimList[3] = { Nres, Nres, Nres };
+  uint32_t outDimList2[4] = { 3, Nres, Nres, Nres };
+  ProgressiveOutput<float> out_den_1 = 
+    ProgressiveOutput<float>::saveArrayProgressive("density.nc", outDimList, 3);
+  ProgressiveOutput<float> out_vel_1 = 
+    ProgressiveOutput<float>::saveArrayProgressive("v3d.nc", outDimList2, 4);
+  ProgressiveOutput<float> out_rad_1 = 
+    ProgressiveOutput<float>::saveArrayProgressive("rad.nc", outDimList, 3);
+ 
+  cout << "Weighing..." << endl;
+  for (int rz = 0; rz < Nres; rz++)
+   {
+      double pz = (rz)*boxsize/Nres-cz;
+
+      cout << rz << " / " << Nres << endl;
+      for (int ry = 0; ry < Nres; ry++)
+	{
+	  double py = (ry)*boxsize/Nres-cy;
+	  for (int rx = 0; rx < Nres; rx++)
+	    {
+	      double px = (rx)*boxsize/Nres-cx;
+	      
+	      MyTree::coords c = { px, py, pz };
+
+	      double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
+	      if (r2 > rLimit2)
+		{
+		  continue;
+		}
+
+	      uint32_t numInCell = bins[rx + Nres*ry + Nres*Nres*rz];
+	      if (numInCell > N_SPH)
+		smooth1.fetchNeighbours(c, numInCell);
+	      else
+		smooth1.fetchNeighbours(c);
+	      smooth1.addGridSite(c);
+	    }
+	}
+    }
+
+  
+  
+
+  cout << "Interpolating..." << endl;
+  for (int rz = 0; rz < Nres; rz++)
+    {
+      double pz = (rz)*boxsize/Nres-cz;
+
+      cout << rz << " / " << Nres << endl;
+      for (int ry = 0; ry < Nres; ry++)
+	{
+	  double py = (ry)*boxsize/Nres-cy;
+	  for (int rx = 0; rx < Nres; rx++)
+	    {
+	      double px = (rx)*boxsize/Nres-cx;
+	      
+	      MyTree::coords c = { px, py, pz };
+
+	      double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
+	      if (r2 > rLimit2)
+		{
+		  out_vel_1.put(0);
+		  out_vel_1.put(0);
+		  out_vel_1.put(0);
+		  out_den_1.put(0);
+		  out_rad_1.put(0);
+		  continue;
+		}
+
+	      uint32_t numInCell = bins[rx + ry*Nres + rz*Nres*Nres];
+	      if (numInCell > N_SPH)
+		smooth1.fetchNeighbours(c, numInCell);
+	      else
+		smooth1.fetchNeighbours(c);
+
+	      float val;
+
+	      out_rad_1.put(smooth1.getSmoothingLen());
+	      val = smooth1.computeSmoothedValue(c, getVelocity<0>);
+	      out_vel_1.put(val);
+	      val = smooth1.computeSmoothedValue(c, getVelocity<1>);
+	      out_vel_1.put(val);
+	      val = smooth1.computeSmoothedValue(c, getVelocity<2>);
+	      out_vel_1.put(val);
+	      val = smooth1.computeSmoothedValue(c, getUnity);
+	      out_den_1.put(val);
+	    }
+	}
+    }
+
+  return 0;
+};

From 34de69986398ff674fadbf8828d9e4e32d4d27ee Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem.lavaux@laposte.net>
Date: Tue, 20 May 2014 09:43:00 +0200
Subject: [PATCH 13/19] Added proper support for boost detection

---
 CMakeLists.txt            |  2 ++
 sample/CMakeLists.txt     | 17 +++++++++++------
 sample/simple3DFilter.cpp |  8 ++++----
 3 files changed, 17 insertions(+), 10 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index e2e45f6..1830463 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -31,6 +31,8 @@ find_library(NETCDF_LIBRARY netcdf)
 find_library(GSL_LIBRARY gsl)
 find_library(GSLCBLAS_LIBRARY gslcblas)
 
+find_package(Boost 1.53)
+
 if (ENABLE_SHARP)
   SET(SHARP_SOURCE ${CMAKE_SOURCE_DIR}/external/sharp)
   SET(DEP_BUILD ${CMAKE_SOURCE_DIR}/external/sharp/auto)
diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt
index e7f60fe..a333aff 100644
--- a/sample/CMakeLists.txt
+++ b/sample/CMakeLists.txt
@@ -1,5 +1,7 @@
 SET(tolink ${GSL_LIBRARIES} CosmoTool ${CosmoTool_LIBS})
-include_directories(${CMAKE_SOURCE_DIR}/src ${FFTW3_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH})
+include_directories(${CMAKE_SOURCE_DIR}/src)
+include_directories(${FFTW3_INCLUDE_DIRS}
+${EIGEN3_INCLUDE_DIRS} ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH})
 
 IF(SHARP_INCLUDE_PATH)
   include_directories(BEFORE ${SHARP_INCLUDE_PATH})
@@ -11,9 +13,6 @@ target_link_libraries(testBQueue ${tolink})
 add_executable(testInterpolate testInterpolate.cpp)
 target_link_libraries(testInterpolate ${tolink})
 
-add_executable(testSmooth testSmooth.cpp)
-target_link_libraries(testSmooth ${tolink})
-
 add_executable(testkd testkd.cpp)
 target_link_libraries(testkd ${tolink})
 
@@ -69,6 +68,12 @@ endif (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
 add_executable(test_cosmopower test_cosmopower.cpp)
 target_link_libraries(test_cosmopower ${tolink})
 
+if (Boost_FOUND)
+  include_directories(${Boost_INCLUDE_DIRS})
 
-add_executable(simple3DFilter simple3DFilter.cpp)
-target_link_libraries(simple3DFilter ${tolink})
+  add_executable(testSmooth testSmooth.cpp)
+  target_link_libraries(testSmooth ${tolink})
+
+  add_executable(simple3DFilter simple3DFilter.cpp)
+  target_link_libraries(simple3DFilter ${tolink})
+endif (Boost_FOUND)
diff --git a/sample/simple3DFilter.cpp b/sample/simple3DFilter.cpp
index 8b91332..b7d1498 100644
--- a/sample/simple3DFilter.cpp
+++ b/sample/simple3DFilter.cpp
@@ -1,8 +1,8 @@
 #include <cassert>
-#include <CosmoTool/yorick.hpp>
-#include <CosmoTool/sphSmooth.hpp>
-#include <CosmoTool/mykdtree.hpp>
-#include <CosmoTool/miniargs.hpp>
+#include "yorick.hpp"
+#include "sphSmooth.hpp"
+#include "mykdtree.hpp"
+#include "miniargs.hpp"
 
 using namespace std;
 using namespace CosmoTool;

From c6668cbd0a6217261494842184dd8b3a844704ba Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem.lavaux@laposte.net>
Date: Tue, 20 May 2014 10:58:07 +0200
Subject: [PATCH 14/19] Fixes to sphSmooth

---
 src/sphSmooth.tcc | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/sphSmooth.tcc b/src/sphSmooth.tcc
index b9c9d2e..1ae5486 100644
--- a/src/sphSmooth.tcc
+++ b/src/sphSmooth.tcc
@@ -47,12 +47,12 @@ SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, uin
   if (requested > maxNgb)
     {
       maxNgb = requested;
-      internal.ngb = new SPHCell *[maxNgb];
-      internal.distances = new CoordType[maxNgb];
+      internal.ngb = boost::shared_ptr<SPHCell*[]>(new SPHCell *[maxNgb]);
+      internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
     }
 
   memcpy(internal.currentCenter, c, sizeof(c));
-  tree->getNearestNeighbours(c, requested, internal.ngb, internal.distances);
+  tree->getNearestNeighbours(c, requested, internal.ngb.get(), internal.distances.get());
 
   internal.currentNgb = 0;
   for (uint32_t i = 0; i < requested && internal.ngb[i] != 0; i++,internal.currentNgb++)

From 11dd968f3677abbbe133091f3930314083f41157 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem.lavaux@laposte.net>
Date: Tue, 20 May 2014 11:19:45 +0200
Subject: [PATCH 15/19] Compile CIC. Imported tool to convert gadget snapshot
 to density field

---
 sample/CMakeLists.txt      |  3 ++
 sample/gadgetToDensity.cpp | 62 ++++++++++++++++++++++++++++++++++++++
 src/CMakeLists.txt         |  1 +
 3 files changed, 66 insertions(+)
 create mode 100644 sample/gadgetToDensity.cpp

diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt
index a333aff..606c8d1 100644
--- a/sample/CMakeLists.txt
+++ b/sample/CMakeLists.txt
@@ -77,3 +77,6 @@ if (Boost_FOUND)
   add_executable(simple3DFilter simple3DFilter.cpp)
   target_link_libraries(simple3DFilter ${tolink})
 endif (Boost_FOUND)
+
+add_executable(gadgetToDensity gadgetToDensity.cpp)
+target_link_libraries(gadgetToDensity ${tolink})
diff --git a/sample/gadgetToDensity.cpp b/sample/gadgetToDensity.cpp
new file mode 100644
index 0000000..d17f263
--- /dev/null
+++ b/sample/gadgetToDensity.cpp
@@ -0,0 +1,62 @@
+#include <cmath>
+#include <iostream>
+#include <cstdlib>
+#include "cic.hpp"
+#include "loadGadget.hpp"
+#include "miniargs.hpp"
+#include "yorick.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+int main(int argc, char **argv)
+{
+  uint32_t res;
+  char *fname;
+  int id;
+
+  MiniArgDesc desc[] = {
+    { "SNAPSHOT", &fname, MINIARG_STRING },
+    { "ID", &id, MINIARG_INT },
+    { "RESOLUTION", &res, MINIARG_INT },
+    { 0, 0, MINIARG_NULL }
+  };
+
+  if (!parseMiniArgs(argc, argv, desc))
+    return 1;
+
+  SimuData *p = loadGadgetMulti(fname, 0, 0);
+  double L0 = p->BoxSize/1000;
+  CICFilter filter(res, L0);
+
+  delete p;
+
+  try {
+  for (int cpuid=0;;cpuid++) {
+    p = loadGadgetMulti(fname, cpuid, NEED_POSITION);
+    for (uint32_t i = 0; i < p->NumPart; i++)
+    {
+      CICParticles a;
+      
+      a.mass = 1.0;
+      a.coords[0] = p->Pos[0][i]/1000;
+      a.coords[1] = p->Pos[1][i]/1000;
+      a.coords[2] = p->Pos[2][i]/1000;
+      filter.putParticles(&a, 1);
+    }
+  
+    delete p;
+  }
+  } catch (const NoSuchFileException& e) {}
+
+  CICType *denField;
+  uint32_t Ntot;
+  filter.getDensityField(denField, Ntot);
+  
+  cout << "L0=" << L0 << endl;
+  cout << "Saving density field" << endl;
+  uint32_t dimList[] = { res, res, res};
+  saveArray("densityField.nc", denField, dimList, 3);
+
+  return 0;
+}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 8901027..611deaa 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -9,6 +9,7 @@ SET(CosmoTool_SRCS
 	miniargs.cpp
 	growthFactor.cpp
 	cosmopower.cpp
+  cic.cpp
 )
 
 IF(FOUND_NETCDF3)

From 47255ea25a4fc5b80463229da8ad553e9f4c6b6c Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Tue, 20 May 2014 16:16:46 +0200
Subject: [PATCH 16/19] HDF5 boost multi array transparent support

---
 sample/CMakeLists.txt |   3 +
 sample/testHDF5.cpp   |  65 ++++++++++++++++
 src/hdf5_array.hpp    | 175 ++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 243 insertions(+)
 create mode 100644 sample/testHDF5.cpp
 create mode 100644 src/hdf5_array.hpp

diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt
index 606c8d1..5180106 100644
--- a/sample/CMakeLists.txt
+++ b/sample/CMakeLists.txt
@@ -34,6 +34,9 @@ target_link_libraries(testPool ${tolink})
 if (HDF5_FOUND)
 	  add_executable(testReadFlash testReadFlash.cpp)
     target_link_libraries(testReadFlash ${tolink})
+    
+    add_executable(testHDF5 testHDF5.cpp)
+    target_link_libraries(testHDF5 ${tolink})
 endif (HDF5_FOUND)
 
 
diff --git a/sample/testHDF5.cpp b/sample/testHDF5.cpp
new file mode 100644
index 0000000..82c3230
--- /dev/null
+++ b/sample/testHDF5.cpp
@@ -0,0 +1,65 @@
+#include <iostream>
+#include <H5Cpp.h>
+#include "hdf5_array.hpp"
+
+using namespace std;
+
+int main()
+{
+  typedef  boost::multi_array<float, 2> array_type;
+  typedef  boost::multi_array<float, 3> array3_type;
+  typedef  boost::multi_array<std::complex<double>, 2> arrayc_type;
+  typedef array_type::index index;
+
+  H5::H5File f("test.h5", H5F_ACC_TRUNC);
+  
+  H5::Group g = f.createGroup("test_group");
+  
+  array_type A(boost::extents[2][3]);
+  array_type B;
+  array3_type C(boost::extents[2][3][4]);
+  arrayc_type D, E;
+  
+  int values = 0;
+  for (index i = 0; i != 2; i++)
+    for (index j = 0; j != 3; j++)
+      A[i][j] = values++;
+      
+  CosmoTool::hdf5_write_array(g, "test_data", A);
+  
+  CosmoTool::hdf5_read_array(g, "test_data", B);
+
+  int verify = 0;
+  for (index i = 0; i != 2; i++)
+    for (index j = 0; j != 3; j++)
+      if (B[i][j] != verify++) {
+        std::cout << "Invalid array content" << endl;
+        abort();
+      }
+  
+  try
+    {
+      CosmoTool::hdf5_read_array(g, "test_data", C);
+      std::cout << "Did not throw InvalidDimensions" << endl;
+      abort();
+    }
+  catch (const CosmoTool::InvalidDimensions&)
+    {}
+  
+  D.resize(boost::extents[2][3]);
+  D = A;
+  
+  CosmoTool::hdf5_write_array(g, "test_data_c", D);
+
+  CosmoTool::hdf5_read_array(g, "test_data_c", E);
+
+  verify = 0;
+  for (index i = 0; i != 2; i++)
+    for (index j = 0; j != 3; j++)
+      if (E[i][j].real() != verify++) {
+        std::cout << "Invalid array content" << endl;
+        abort();
+      }
+
+  return 0;
+}
diff --git a/src/hdf5_array.hpp b/src/hdf5_array.hpp
new file mode 100644
index 0000000..46673ba
--- /dev/null
+++ b/src/hdf5_array.hpp
@@ -0,0 +1,175 @@
+#ifndef __COSMO_HDF5_ARRAY_HPP
+#define __COSMO_HDF5_ARRAY_HPP
+
+
+#include <string>
+#include <stdint.h>
+#include <boost/static_assert.hpp>
+#include <boost/utility.hpp>
+#include <boost/mpl/assert.hpp>
+#include <boost/multi_array.hpp>
+#include <boost/type_traits/is_same.hpp> 
+#include <boost/type_traits/is_floating_point.hpp> 
+#include <vector>
+#include <H5Cpp.h>
+
+namespace CosmoTool {
+
+  //!_______________________________________________________________________________________
+  //!     
+  //!     map types to HDF5 types
+  //!         
+  //!     
+  //!     Leo Goodstadt (04 March 2013), improved with enable_if by Guilhem Lavaux (May 2014)
+  //!_______________________________________________________________________________________ 
+
+  template<typename T, class Enable = void> struct get_hdf5_data_type
+  {
+    static H5::PredType type()  
+      {   
+          BOOST_MPL_ASSERT_MSG(0, Unknown_HDF5_data_type, ()); 
+          return H5::PredType::NATIVE_DOUBLE; 
+      }
+  };
+  
+  #define HDF5_TYPE(tl, thdf5) \
+  template<typename T> struct get_hdf5_data_type<T, typename boost::enable_if<boost::is_same<T, tl> >::type>  \
+    { static H5::PredType type() { return H5::PredType::thdf5; }; }
+  
+  HDF5_TYPE(char              , NATIVE_CHAR);
+  HDF5_TYPE(long long         , NATIVE_LLONG);
+  HDF5_TYPE(unsigned long long, NATIVE_ULLONG);
+  HDF5_TYPE(int8_t            , NATIVE_INT8);
+  HDF5_TYPE(uint8_t           , NATIVE_UINT8);
+  HDF5_TYPE(int16_t           , NATIVE_INT16);
+  HDF5_TYPE(uint16_t          , NATIVE_UINT16);
+  HDF5_TYPE(int32_t           , NATIVE_INT32);
+  HDF5_TYPE(uint32_t          , NATIVE_UINT32);
+  HDF5_TYPE(int64_t           , NATIVE_INT64);
+  HDF5_TYPE(uint64_t          , NATIVE_UINT64);
+  HDF5_TYPE(float             , NATIVE_FLOAT);
+  HDF5_TYPE(double            , NATIVE_DOUBLE);
+  HDF5_TYPE(long double       , NATIVE_LDOUBLE);
+  
+  #undef HDF5_TYPE
+
+//!_______________________________________________________________________________________
+//!     
+//!     write_hdf5 multi_array
+//!         
+//!     \author leo Goodstadt (04 March 2013)
+//!     
+//!_______________________________________________________________________________________
+  template<typename T, std::size_t DIMENSIONS, typename hdf5_data_type>
+  void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, 
+                        const boost::multi_array<T, DIMENSIONS>& data, 
+                        const hdf5_data_type& datatype)
+  {
+      std::vector<hsize_t> dimensions(data.shape(), data.shape() + DIMENSIONS);
+      H5::DataSpace dataspace(DIMENSIONS, dimensions.data());
+
+      H5::DataSet dataset = fg.createDataSet(data_set_name, datatype, dataspace);
+
+      dataset.write(data.data(), datatype);
+  }
+  
+  /* HDF5 complex type */
+  template<typename T>
+  class hdf5_ComplexType
+  {
+  public:
+    H5::CompType type;
+    
+    hdf5_ComplexType()
+      : type(sizeof(std::complex<T>))
+    {
+      get_hdf5_data_type<T> hdf_data_type;
+      type.insertMember("r", 0, hdf_data_type.type());
+      type.insertMember("i", sizeof(T), hdf_data_type.type());
+      type.pack();
+    }
+    
+    static const hdf5_ComplexType<T> *ctype()
+    {
+      static hdf5_ComplexType<T> singleton;
+      
+      return &singleton;
+    }
+  };
+
+  template<typename T, std::size_t DIMENSIONS>
+  void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const boost::multi_array<T, DIMENSIONS>& data )
+  {
+      get_hdf5_data_type<T> hdf_data_type;
+      hdf5_write_array(fg, data_set_name, data, hdf_data_type.type());
+  }
+
+  template<typename T, std::size_t DIMENSIONS>
+  void hdf5_write_array(H5::CommonFG& fg, const std::string& data_set_name, const boost::multi_array<std::complex<T>, DIMENSIONS>& data )
+  {
+      hdf5_write_array(fg, data_set_name, data, hdf5_ComplexType<T>::ctype()->type);
+  }
+
+  
+  // HDF5 array reader
+  //
+  // Author Guilhem Lavaux (May 2014)
+  
+  class InvalidDimensions: virtual std::exception {
+  };
+  
+  template<std::size_t r>
+  struct hdf5_extent_gen {
+    typedef typename boost::detail::multi_array::extent_gen<r> type;
+    
+    static inline type build(hsize_t *d)
+    {
+      return (hdf5_extent_gen<r-1>::build(d))[d[r-1]];
+    }
+  };
+
+  template<>
+  struct hdf5_extent_gen<0> {    
+    static inline boost::multi_array_types::extent_gen build(hsize_t *d)
+    {
+      return boost::extents;
+    }
+  };
+
+  template<typename T, std::size_t DIMENSIONS, typename hdf5_data_type>
+  void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, 
+                       boost::multi_array<T, DIMENSIONS>& data, 
+                       const hdf5_data_type& datatype)
+  {
+      H5::DataSet dataset = fg.openDataSet(data_set_name);
+      H5::DataSpace dataspace = dataset.getSpace();
+      std::vector<hsize_t> dimensions(DIMENSIONS);
+      
+      if (dataspace.getSimpleExtentNdims() != DIMENSIONS)
+        {
+          throw InvalidDimensions();
+        }
+            
+      dataspace.getSimpleExtentDims(dimensions.data());
+      data.resize(hdf5_extent_gen<DIMENSIONS>::build(dimensions.data()));
+      dataset.read(data.data(), datatype);
+  }
+
+  template<typename T, std::size_t DIMENSIONS>
+  void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, boost::multi_array<T, DIMENSIONS>& data )
+  {
+      get_hdf5_data_type<T> hdf_data_type;
+      hdf5_read_array(fg, data_set_name, data, hdf_data_type.type());
+  }
+
+  template<typename T, std::size_t DIMENSIONS>
+  void hdf5_read_array(H5::CommonFG& fg, const std::string& data_set_name, boost::multi_array<std::complex<T>, DIMENSIONS>& data )
+  {
+      hdf5_read_array(fg, data_set_name, data, hdf5_ComplexType<T>::ctype()->type);
+  }
+
+};
+
+#endif
+
+

From 468f5222b8f04aef6d51ec0c95ef32aba56e8ba8 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Tue, 20 May 2014 16:27:41 +0200
Subject: [PATCH 17/19] New tool to convert a gadget to a HDF5 file

---
 sample/CMakeLists.txt    |  4 +++
 sample/gadgetToArray.cpp | 54 ++++++++++++++++++++++++++++++++++++++++
 2 files changed, 58 insertions(+)
 create mode 100644 sample/gadgetToArray.cpp

diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt
index 5180106..47a74b5 100644
--- a/sample/CMakeLists.txt
+++ b/sample/CMakeLists.txt
@@ -37,6 +37,10 @@ if (HDF5_FOUND)
     
     add_executable(testHDF5 testHDF5.cpp)
     target_link_libraries(testHDF5 ${tolink})
+
+    add_executable(gadgetToArray gadgetToArray.cpp)
+    target_link_libraries(gadgetToArray ${tolink})
+
 endif (HDF5_FOUND)
 
 
diff --git a/sample/gadgetToArray.cpp b/sample/gadgetToArray.cpp
new file mode 100644
index 0000000..eb01e9a
--- /dev/null
+++ b/sample/gadgetToArray.cpp
@@ -0,0 +1,54 @@
+#include <cmath>
+#include <iostream>
+#include <cstdlib>
+#include "cic.hpp"
+#include "loadGadget.hpp"
+#include "miniargs.hpp"
+#include <H5Cpp.h>
+#include "hdf5_array.hpp"
+
+using namespace CosmoTool;
+using namespace std;
+
+int main(int argc, char **argv)
+{
+  typedef boost::multi_array<float, 2> array_type;
+  uint32_t res;
+  char *fname;
+  int id;
+
+  MiniArgDesc desc[] = {
+    { "SNAPSHOT", &fname, MINIARG_STRING },
+    { 0, 0, MINIARG_NULL }
+  };
+
+  if (!parseMiniArgs(argc, argv, desc))
+    return 1;
+
+  H5::H5File f("density.h5", H5F_ACC_TRUNC);
+  
+
+  SimuData *p = loadGadgetMulti(fname, 0, 0);
+  double L0 = p->BoxSize/1000;
+  array_type parts(boost::extents[p->TotalNumPart][3]);
+  uint64_t q = 0;
+  
+  try {
+    for (int cpuid=0;;cpuid++) {
+      p = loadGadgetMulti(fname, cpuid, NEED_POSITION);
+      for (uint32_t i = 0; i < p->NumPart; i++)
+        {
+          parts[q][0] = p->Pos[0][i]/1000;
+          parts[q][1] = p->Pos[1][i]/1000;
+          parts[q][2] = p->Pos[2][i]/1000;
+          q++;
+        }
+
+      delete p;
+    }
+  } catch (const NoSuchFileException& e) {}
+
+  hdf5_write_array(f, "particles", parts);
+
+  return 0;
+}

From c88599aa594612f5f641ddc5b5c76f41a6e43e01 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Tue, 20 May 2014 16:53:51 +0200
Subject: [PATCH 18/19] Export also velocities

---
 sample/gadgetToArray.cpp | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/sample/gadgetToArray.cpp b/sample/gadgetToArray.cpp
index eb01e9a..7d4669c 100644
--- a/sample/gadgetToArray.cpp
+++ b/sample/gadgetToArray.cpp
@@ -30,17 +30,20 @@ int main(int argc, char **argv)
 
   SimuData *p = loadGadgetMulti(fname, 0, 0);
   double L0 = p->BoxSize/1000;
-  array_type parts(boost::extents[p->TotalNumPart][3]);
+  array_type parts(boost::extents[p->TotalNumPart][6]);
   uint64_t q = 0;
   
   try {
     for (int cpuid=0;;cpuid++) {
-      p = loadGadgetMulti(fname, cpuid, NEED_POSITION);
+      p = loadGadgetMulti(fname, cpuid, NEED_POSITION|NEED_VELOCITY);
       for (uint32_t i = 0; i < p->NumPart; i++)
         {
           parts[q][0] = p->Pos[0][i]/1000;
           parts[q][1] = p->Pos[1][i]/1000;
           parts[q][2] = p->Pos[2][i]/1000;
+          parts[q][3] = p->Vel[0][i];
+          parts[q][4] = p->Vel[1][i];
+          parts[q][5] = p->Vel[2][i];
           q++;
         }
 

From 113a14e00ca1499bf7d6980ed5faab6a8a5fa9c5 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Tue, 20 May 2014 17:16:54 +0200
Subject: [PATCH 19/19] Parallelized and HDF5ized simple3DFilter

---
 sample/simple3DFilter.cpp | 199 ++++++++++++++++++++------------------
 1 file changed, 105 insertions(+), 94 deletions(-)

diff --git a/sample/simple3DFilter.cpp b/sample/simple3DFilter.cpp
index b7d1498..1932e97 100644
--- a/sample/simple3DFilter.cpp
+++ b/sample/simple3DFilter.cpp
@@ -3,6 +3,9 @@
 #include "sphSmooth.hpp"
 #include "mykdtree.hpp"
 #include "miniargs.hpp"
+#include <H5Cpp.h>
+#include "hdf5_array.hpp"
+#include <iostream>
 
 using namespace std;
 using namespace CosmoTool;
@@ -30,6 +33,10 @@ typedef MyTree::Cell MyCell;
 
 int main(int argc, char **argv)
 {
+  typedef boost::multi_array<float, 2> array_type;
+  typedef boost::multi_array<float, 3> array3_type;
+  typedef boost::multi_array<float, 4> array4_type;
+
   char *fname1, *fname2;
   double rLimit, boxsize, rLimit2, cx, cy, cz;
   int Nres;
@@ -48,35 +55,34 @@ int main(int argc, char **argv)
   if (!parseMiniArgs(argc, argv, args))
     return 1;
 
-  float *v1_data;
-  uint32_t *dimList;
-  uint32_t rank;
+  H5::H5File in_f(fname1, 0);
+  H5::H5File out_f("fields.h5", H5F_ACC_TRUNC);
+  array_type v1_data;
   uint32_t N1_points, N2_points;
-  int *bins = new int[Nres*Nres*Nres];
+  
+  array3_type bins(boost::extents[Nres][Nres][Nres]);
 
   rLimit2 = rLimit*rLimit;
 
-  loadArray(fname1, v1_data, dimList, rank);
-  assert(rank == 2);
-  assert(dimList[1] == 6);
+  hdf5_read_array(in_f, "particles", v1_data);
+  assert(v1_data.shape()[1] == 6);
 
-  N1_points = dimList[0];
-  delete[] dimList;
+  N1_points = v1_data.shape()[0];
   
   cout << "Got " << N1_points << " in the first file." << endl;
 
   MyCell *allCells_1 = new MyCell[N1_points];
   
   for (long i = 0; i < Nres*Nres*Nres; i++)
-    bins[i] = 0;
+    bins.data()[i] = 0;
 
   cout << "Shuffling data in cells..." << endl;
   for (int i = 0 ; i < N1_points; i++)
     {
       for (int j = 0; j < 3; j++)
-	allCells_1[i].coord[j] = v1_data[i*6 + j];
+        allCells_1[i].coord[j] = v1_data[i][j];
       for (int k = 0; k < 3; k++)
-	allCells_1[i].val.pValue.v[k] = v1_data[i*6 + 3 + k];
+        allCells_1[i].val.pValue.v[k] = v1_data[i][3+k];
       allCells_1[i].active = true;
       allCells_1[i].val.weight = 0.0;
 
@@ -85,109 +91,114 @@ int main(int argc, char **argv)
       long rz = floor((allCells_1[i].coord[2]+cz)*Nres/boxsize+0.5);
       
       if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres)
-	  continue;
+        continue;
 	  
-      bins[rx + ry*Nres + rz*Nres*Nres]++;
+      bins[rx][ry][rz]++;
     }
-  delete[] v1_data;
-   uint32_t dims[3] = { Nres, Nres, Nres } ;
-   saveArray("num_in_cell.nc", bins, dims, 3);
+  v1_data.resize(boost::extents[1][1]);
+  
+  hdf5_write_array(out_f, "num_in_cell", bins);
 
   cout << "Building trees..." << endl;
   MyTree tree1(allCells_1, N1_points);
 
   cout << "Creating smoothing filter..." << endl;
-  MySmooth smooth1(&tree1, N_SPH);
 
-  uint32_t outDimList[3] = { Nres, Nres, Nres };
-  uint32_t outDimList2[4] = { 3, Nres, Nres, Nres };
-  ProgressiveOutput<float> out_den_1 = 
-    ProgressiveOutput<float>::saveArrayProgressive("density.nc", outDimList, 3);
-  ProgressiveOutput<float> out_vel_1 = 
-    ProgressiveOutput<float>::saveArrayProgressive("v3d.nc", outDimList2, 4);
-  ProgressiveOutput<float> out_rad_1 = 
-    ProgressiveOutput<float>::saveArrayProgressive("rad.nc", outDimList, 3);
+  array3_type out_den_1(boost::extents[Nres][Nres][Nres]);
+  array4_type out_v3d_1(boost::extents[Nres][Nres][Nres][3]);
+  array3_type out_rad_1(boost::extents[Nres][Nres][Nres]);
  
   cout << "Weighing..." << endl;
-  for (int rz = 0; rz < Nres; rz++)
-   {
-      double pz = (rz)*boxsize/Nres-cz;
 
-      cout << rz << " / " << Nres << endl;
-      for (int ry = 0; ry < Nres; ry++)
-	{
-	  double py = (ry)*boxsize/Nres-cy;
-	  for (int rx = 0; rx < Nres; rx++)
-	    {
-	      double px = (rx)*boxsize/Nres-cx;
-	      
-	      MyTree::coords c = { px, py, pz };
+#pragma omp parallel
+  {
+    MySmooth smooth1(&tree1, N_SPH);
+    
+#pragma omp for schedule(dynamic) 
+    for (int rz = 0; rz < Nres; rz++)
+     {
+        double pz = (rz)*boxsize/Nres-cz;
 
-	      double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
-	      if (r2 > rLimit2)
-		{
-		  continue;
-		}
+        cout << rz << " / " << Nres << endl;
+        for (int ry = 0; ry < Nres; ry++)
+          {
+            double py = (ry)*boxsize/Nres-cy;
+            for (int rx = 0; rx < Nres; rx++)
+              {
+                double px = (rx)*boxsize/Nres-cx;
+                
+                MyTree::coords c = { px, py, pz };
 
-	      uint32_t numInCell = bins[rx + Nres*ry + Nres*Nres*rz];
-	      if (numInCell > N_SPH)
-		smooth1.fetchNeighbours(c, numInCell);
-	      else
-		smooth1.fetchNeighbours(c);
-	      smooth1.addGridSite(c);
-	    }
-	}
-    }
+                double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
+                if (r2 > rLimit2)
+                  {
+                    continue;
+                  }
 
+                uint32_t numInCell = bins[rx][ry][rz];
+                if (numInCell > N_SPH)
+                  smooth1.fetchNeighbours(c, numInCell);
+                else
+                  smooth1.fetchNeighbours(c);
+#pragma omp critical
+                smooth1.addGridSite(c);
+              }
+          }
+      }  
+  }
   
-  
-
   cout << "Interpolating..." << endl;
-  for (int rz = 0; rz < Nres; rz++)
-    {
-      double pz = (rz)*boxsize/Nres-cz;
+#pragma omp parallel
+  {
+    MySmooth smooth1(&tree1, N_SPH);
+    
+#pragma omp for schedule(dynamic) 
+    for (int rz = 0; rz < Nres; rz++)
+      {
+        double pz = (rz)*boxsize/Nres-cz;
 
-      cout << rz << " / " << Nres << endl;
-      for (int ry = 0; ry < Nres; ry++)
-	{
-	  double py = (ry)*boxsize/Nres-cy;
-	  for (int rx = 0; rx < Nres; rx++)
-	    {
-	      double px = (rx)*boxsize/Nres-cx;
-	      
-	      MyTree::coords c = { px, py, pz };
+        cout << rz << " / " << Nres << endl;
+        for (int ry = 0; ry < Nres; ry++)
+          {
+            double py = (ry)*boxsize/Nres-cy;
+            for (int rx = 0; rx < Nres; rx++)
+              {
+                double px = (rx)*boxsize/Nres-cx;
+                
+                MyTree::coords c = { px, py, pz };
 
-	      double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
-	      if (r2 > rLimit2)
-		{
-		  out_vel_1.put(0);
-		  out_vel_1.put(0);
-		  out_vel_1.put(0);
-		  out_den_1.put(0);
-		  out_rad_1.put(0);
-		  continue;
-		}
+                double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
+                if (r2 > rLimit2)
+                  {
+                    out_v3d_1[rx][ry][rz][0] = 0;
+                    out_v3d_1[rx][ry][rz][1] = 0;
+                    out_v3d_1[rx][ry][rz][2] = 0;
+                    out_den_1[rx][ry][rz] = 0;
+                    out_rad_1[rx][ry][rz] = -1;
+                    continue;
+                  }
 
-	      uint32_t numInCell = bins[rx + ry*Nres + rz*Nres*Nres];
-	      if (numInCell > N_SPH)
-		smooth1.fetchNeighbours(c, numInCell);
-	      else
-		smooth1.fetchNeighbours(c);
+                uint32_t numInCell = bins[rx][ry][rz];
+                if (numInCell > N_SPH)
+                  smooth1.fetchNeighbours(c, numInCell);
+                else
+                  smooth1.fetchNeighbours(c);
 
-	      float val;
+                float val;
 
-	      out_rad_1.put(smooth1.getSmoothingLen());
-	      val = smooth1.computeSmoothedValue(c, getVelocity<0>);
-	      out_vel_1.put(val);
-	      val = smooth1.computeSmoothedValue(c, getVelocity<1>);
-	      out_vel_1.put(val);
-	      val = smooth1.computeSmoothedValue(c, getVelocity<2>);
-	      out_vel_1.put(val);
-	      val = smooth1.computeSmoothedValue(c, getUnity);
-	      out_den_1.put(val);
-	    }
-	}
-    }
+                out_rad_1[rx][ry][rz] = smooth1.getSmoothingLen();
+                out_v3d_1[rx][ry][rz][0] = smooth1.computeSmoothedValue(c, getVelocity<0>);
+                out_v3d_1[rx][ry][rz][1] = smooth1.computeSmoothedValue(c, getVelocity<1>);
+                out_v3d_1[rx][ry][rz][2] = smooth1.computeSmoothedValue(c, getVelocity<2>);
+                out_den_1[rx][ry][rz] = smooth1.computeSmoothedValue(c, getUnity);
+              }
+          }
+      }
+  }
+      
+  hdf5_write_array(out_f, "radii", out_rad_1);
+  hdf5_write_array(out_f, "velocity", out_v3d_1);
+  hdf5_write_array(out_f, "density", out_den_1);
 
   return 0;
 };