Merge branch 'master' of bitbucket.org:glavaux/cosmotool

CMakeLists.txt

@@ -55,6 +55,7 @@ include(FindHDF5)
 include(FindPkgConfig)
 
 pkg_check_modules(FFTW3 fftw3>=3.3)
+pkg_check_modules(FFTW3F fftw3f>=3.3)
 pkg_check_modules(EIGEN3 eigen3)
 
 include(FindPackageHandleStandardArgs)

build_tools/gather_sources.py

@@ -0,0 +1,161 @@
+#+
+# This is CosmoTool (./build_tools/gather_sources.py) -- Copyright (C) Guilhem Lavaux (2007-2013)
+#
+# guilhem.lavaux@gmail.com
+#
+# This software is a computer program whose purpose is to provide a toolbox for cosmological
+# data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+#
+# This software is governed by the CeCILL license under French law and
+# abiding by the rules of distribution of free software.  You can  use, 
+# modify and/ or redistribute the software under the terms of the CeCILL
+# license as circulated by CEA, CNRS and INRIA at the following URL
+# "http://www.cecill.info". 
+#
+# As a counterpart to the access to the source code and  rights to copy,
+# modify and redistribute granted by the license, users are provided only
+# with a limited warranty  and the software's author,  the holder of the
+# economic rights,  and the successive licensors  have only  limited
+# liability. 
+#
+# In this respect, the user's attention is drawn to the risks associated
+# with loading,  using,  modifying and/or developing or reproducing the
+# software by the user in light of its specific status of free software,
+# that may mean  that it is complicated to manipulate,  and  that  also
+# therefore means  that it is reserved for developers  and  experienced
+# professionals having in-depth computer knowledge. Users are therefore
+# encouraged to load and test the software's suitability as regards their
+# requirements in conditions enabling the security of their systems and/or 
+# data to be ensured and,  more generally, to use and operate it in the 
+# same conditions as regards security. 
+#
+# The fact that you are presently reading this means that you have had
+# knowledge of the CeCILL license and that you accept its terms.
+#+
+import shutil
+import tempfile
+import re
+from git import Repo,Tree,Blob
+
+def apply_license(license, relimit, filename):
+  header = re.sub(r'@FILENAME@', filename, license)
+
+  f = file(filename)
+  lines = f.read()
+  f.close()
+
+  lines = re.sub(relimit, '', lines)
+  lines = header + lines
+
+  with tempfile.NamedTemporaryFile(delete=False) as tmp_sources:
+    tmp_sources.write(lines)
+
+  shutil.move(tmp_sources.name, filename)
+
+
+def apply_python_license(filename):
+  license="""#+
+# This is CosmoTool (@FILENAME@) -- Copyright (C) Guilhem Lavaux (2007-2013)
+#
+# guilhem.lavaux@gmail.com
+#
+# This software is a computer program whose purpose is to provide a toolbox for cosmological
+# data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+#
+# This software is governed by the CeCILL license under French law and
+# abiding by the rules of distribution of free software.  You can  use, 
+# modify and/ or redistribute the software under the terms of the CeCILL
+# license as circulated by CEA, CNRS and INRIA at the following URL
+# "http://www.cecill.info". 
+#
+# As a counterpart to the access to the source code and  rights to copy,
+# modify and redistribute granted by the license, users are provided only
+# with a limited warranty  and the software's author,  the holder of the
+# economic rights,  and the successive licensors  have only  limited
+# liability. 
+#
+# In this respect, the user's attention is drawn to the risks associated
+# with loading,  using,  modifying and/or developing or reproducing the
+# software by the user in light of its specific status of free software,
+# that may mean  that it is complicated to manipulate,  and  that  also
+# therefore means  that it is reserved for developers  and  experienced
+# professionals having in-depth computer knowledge. Users are therefore
+# encouraged to load and test the software's suitability as regards their
+# requirements in conditions enabling the security of their systems and/or 
+# data to be ensured and,  more generally, to use and operate it in the 
+# same conditions as regards security. 
+#
+# The fact that you are presently reading this means that you have had
+# knowledge of the CeCILL license and that you accept its terms.
+#+
+"""
+   
+  print("Shell/Python file: %s" % filename)
+  relimit=r'^#\+\n(#.*\n)*#\+\n'
+  apply_license(license, relimit, filename)
+
+
+def apply_cpp_license(filename):
+  license="""/*+
+This is CosmoTool (@FILENAME@) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+"""
+  relimit = r'^(?s)/\*\+.*\+\*/'
+  print("C++ file: %s" % filename)
+  apply_license(license, relimit, filename)
+  
+
+def analyze_tree(prefix, t):
+  for ename,entry in t.items():
+    if ename == 'external' or ename == 'zobov':
+      continue
+    if type(entry) == Tree:
+      analyze_tree(prefix + "/" + ename, entry)
+    elif type(entry) == Blob:
+      if ename == './src/hdf5_flash.h' or ename == './src/h5_readFlash.cpp' or ename == './src/h5_readFlash.hpp':
+        continue
+
+      if re.match(".*\.(sh|py|pyx)$",ename) != None:
+        fname=prefix+"/"+ename
+        apply_python_license(fname)
+      if re.match('.*\.(cpp|hpp|h)$', ename) != None:
+        fname=prefix+"/"+ename
+        apply_cpp_license(fname)
+
+
+if __name__=="__main__":
+  repo = Repo(".")
+  assert repo.bare == False
+  t = repo.tree()
+  analyze_tree(".", t)

sample/CMakeLists.txt

@@ -20,6 +20,9 @@ target_link_libraries(testkd ${tolink})
 add_executable(testkd2 testkd2.cpp)
 target_link_libraries(testkd2 ${tolink})
 
+add_executable(testkd3 testkd3.cpp)
+target_link_libraries(testkd3 ${tolink})
+
 add_executable(testDelaunay testDelaunay.cpp)
 target_link_libraries(testDelaunay ${tolink})
 
@@ -44,10 +47,10 @@ target_link_libraries(testAlgo ${tolink})
 add_executable(testBSP testBSP.cpp)
 target_link_libraries(testBSP ${tolink})
 
-if (FFTW3_FOUND AND EIGEN3_FOUND)
+if (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND)
   add_executable(test_fft_calls test_fft_calls.cpp)
-  target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES})
+  target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES})
-endif (FFTW3_FOUND AND EIGEN3_FOUND)
+endif (FFTW3_FOUND AND FFTW3F_FOUND AND  EIGEN3_FOUND)
 
 if (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
   include_directories(${SHARP_INCLUDE_PATH})
@@ -55,3 +58,7 @@ if (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
   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})
 endif (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
+
+add_executable(test_cosmopower test_cosmopower.cpp)
+target_link_libraries(test_cosmopower ${tolink})
+

sample/testAlgo.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/testAlgo.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <iostream>
 #include <stdio.h>
 #include "algo.hpp"

sample/testBQueue.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/testBQueue.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <iostream>
 #include "bqueue.hpp"
 

sample/testBSP.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/testBSP.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #if 0
 #include "bsp_simple.hpp"
 

sample/testDelaunay.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/testDelaunay.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <fstream>
 #include <iostream>
 #include "dinterpolate.hpp"

sample/testEskow.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/testEskow.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <fstream>
 #include <cstring>
 #include <iostream>

sample/testInterpolate.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/testInterpolate.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <iostream>
 #include "interpolate3d.hpp"
 

sample/testNewton.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/testNewton.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <iostream>
 #include <cmath>
 

sample/testPool.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/testPool.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include "pool.hpp"
 
 using namespace CosmoTool;

sample/testReadFlash.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/testReadFlash.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <iostream>
 #include <assert.h>
 #include <stdio.h>

sample/testSmooth.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/testSmooth.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <iostream>
 #include "sphSmooth.hpp"
 #include "yorick.hpp"

sample/test_cosmopower.cpp

@@ -0,0 +1,63 @@
+/*+
+This is CosmoTool (./sample/test_cosmopower.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+#include <cmath>
+#include <iostream>
+#include "cosmopower.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+int main()
+{
+   CosmoPower cp;
+   CosmoPower::CosmoFunction f[] = { CosmoPower::POWER_EFSTATHIOU, CosmoPower::HU_WIGGLES, CosmoPower::HU_BARYON, CosmoPower::POWER_SUGIYAMA };
+   int num_F = sizeof(f)/sizeof(f[0]);
+
+   cp.setFunction(f[0]);
+   cp.normalize();
+   for (int ik = 0; ik < 100; ik++)
+    {
+      double k = pow(10.0, 4*ik/100.-3);
+      
+      cout << k << " ";
+      for (int q = 0; q < num_F; q++)
+        {
+          cp.setFunction(f[q]);
+          cout << cp.power(k) << " ";
+        } 
+      cout << endl;
+    }
+   return 0;
+}

sample/test_fft_calls.cpp

@@ -1,20 +1,93 @@
+/*+
+This is CosmoTool (./sample/test_fft_calls.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#include <boost/format.hpp>
+#include "yorick.hpp"
+#include <gsl/gsl_rng.h>
 #include <iostream>
+#include <cmath>
 #include "fourier/euclidian.hpp"
 
 using namespace CosmoTool;
+using boost::format;
+using boost::str;
 using namespace std;
 
-int main()
+double spectrum_generator(double k)
 {
-  EuclidianFourierTransform_2d<double> dft(128,128,1.0,1.0);
+  if (k==0)
-  double volume = 128*128;
+    return 0;
+  return 1/(0.01+pow(k, 3.0));
+}
+
+template<typename T>
+void test_2d(int Nx, int Ny)
+{
+  EuclidianFourierTransform_2d<T> dft(Nx,Ny,1.0,1.0);
+  EuclidianSpectrum_1D<T> spectrum(spectrum_generator);
+  double volume = Nx*Ny;
+  gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
 
   dft.realSpace().eigen().setRandom();
+  dft.fourierSpace().scale(std::complex<T>(0,0));
   dft.analysis();
+  cout << format("Testing (%d,%d)") % Nx % Ny << endl;
   cout << "Map dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl;
   cout << "Fourier dot-product = " << dft.fourierSpace().dot_product(dft.fourierSpace()).real()*volume << endl;
   dft.synthesis();
   cout << "Resynthesis dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl;
-  return 0;
 
+  dft.realSpace().scale(2.0);
+  dft.fourierSpace().scale(0.2);
+
+  spectrum.newRandomFourier(rng, dft.fourierSpace()); 
+  dft.synthesis();
+
+  uint32_t dims[2] = { Ny, Nx };
+  CosmoTool::saveArray(str(format("generated_map_%d_%d.nc") %Nx % Ny) , dft.realSpace().data(), dims, 2);
+
+  gsl_rng_free(rng);
+}
+
+int main(int argc, char **argv)
+{
+  test_2d<double>(128,128);
+  test_2d<double>(131,128);  
+  test_2d<double>(130,128);  
+  test_2d<float>(128,128);
+  test_2d<float>(128,131);  
+  test_2d<float>(128,130);  
+  return 0;
 }

sample/test_healpix_calls.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/test_healpix_calls.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <iostream>
 #include "fourier/healpix.hpp"
 
@@ -6,15 +41,20 @@ using namespace std;
 
 int main()
 {
-  HealpixFourierTransform<double> dft(128,3*128,3*128, 40);
+  HealpixFourierTransform<double> dft(128,2*128,2*128, 40);
+  HealpixUtilityFunction<double> utils;
   long Npix = dft.realSpace().size();
 
   dft.realSpace().eigen().setRandom();
   dft.analysis();
   cout << "Map dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl;
   cout << "Fourier dot-product = " << dft.fourierSpace().dot_product(dft.fourierSpace()).real()*Npix/(4*M_PI) << endl;
+  cout << "Spectrum analysis" << endl;
+  HealpixUtilityFunction<double>::Spectrum_ptr s = utils.estimateSpectrumFromMap(dft.fourierSpace());
+
   dft.synthesis();
   cout << "Resynthesis dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl;
+
   return 0;
 
 }

sample/testkd.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/testkd.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #define __KD_TREE_SAVE_ON_DISK
 #include <ctime>
 #include <cassert>

sample/testkd2.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./sample/testkd2.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <ctime>
 #include <cassert>
 #include <cstdlib>

sample/testkd3.cpp

@@ -0,0 +1,187 @@
+/*+
+This is CosmoTool (./sample/testkd3.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#include <ctime>
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include <fstream>
+#define __KD_TREE_NUMNODES
+#include "mykdtree.hpp"
+#include "kdtree_splitters.hpp"
+#include <limits>
+
+#define NTRY 30
+#define NGB 24
+#define ND 3
+
+using namespace std;
+using namespace CosmoTool;
+
+typedef KDTree<ND,char,ComputePrecision,KD_homogeneous_cell_splitter<ND, char> > MyTree;
+//typedef KDTree<ND,char,ComputePrecision > MyTree;
+typedef KDCell<ND,char> MyCell;
+
+static double periodic(double a)
+{
+  while (a < -0.5)
+    a+=1;
+  while (a > 0.5)
+    a -= 1;
+  return a;
+}
+
+MyCell *findNearest(MyTree::coords& xc, MyCell *cells, uint32_t Ncells)
+{
+    MyCell *near2 = 0;
+    double R2 = INFINITY;
+    for (int i = 0; i < Ncells; i++)
+      {
+	double d2 = 0;
+	for (int j = 0; j < ND; j++)
+	  {
+	    double delta = periodic(xc[j]-cells[i].coord[j]);
+	    d2 += delta*delta;
+	  }
+	if (d2 < R2)
+	  {
+	    near2 = &cells[i];
+	    R2 = d2;
+	  }
+      }
+    return near2;
+}
+
+int main()
+{
+  uint32_t Ncells = 3000000;
+  MyCell *cells = new MyCell[Ncells];
+  
+  for (int i = 0; i < Ncells; i++)
+    {
+      cells[i].active = true;
+      for (int l = 0; l < ND; l++)
+	cells[i].coord[l] = drand48();
+    }
+  
+  // Check timing
+  clock_t startTimer = clock();
+  MyTree tree(cells, Ncells);
+  clock_t endTimer = clock();
+
+  tree.setPeriodic(true, 1.0);
+
+  clock_t delta = endTimer-startTimer;
+  double myTime = delta*1.0/CLOCKS_PER_SEC * 1.0;
+
+  cout << "KDTree build = " << myTime << " s" << endl;
+
+  MyTree::coords *xc = new MyTree::coords[NTRY];
+
+  cout << "Generating seeds..." << endl;
+  for (int k = 0; k < NTRY; k++)
+    {
+      for (int l = 0; l < ND; l++)
+	xc[k][l] = drand48();
+    }
+
+  // Check consistency
+  cout << "Check consistency..." << endl;
+  MyCell **ngb = new MyCell *[NGB];
+  double *distances = new double[NGB];
+  
+  ofstream fngb("nearest.txt");
+  for (int k = 0; k < NTRY; k++) {
+    cout << "Seed = " << xc[k][0] << " " << xc[k][1] << " " << xc[k][2] << endl;
+    tree.getNearestNeighbours(xc[k], NGB, ngb, distances);
+    int last = -1;
+
+    for (uint32_t i = 0; i < NGB; i++)
+      {
+	if (ngb[i] == 0)
+	  continue;
+
+	last = i;
+
+	double d2 = 0;
+	for (int l = 0; l < 3; l++)
+	  d2 += ({double delta = periodic(xc[k][l] - ngb[i]->coord[l]); delta*delta;});
+	fngb << ngb[i]->coord[0] << " " << ngb[i]->coord[1] << " " << ngb[i]->coord[2] << " " << sqrt(d2) << " " << sqrt(distances[i]) << endl;
+      }
+    fngb << endl << endl;
+    double farther_dist = distances[last];
+    for (uint32_t i = 0; i < Ncells; i++)
+      {
+         bool found = false;
+         // If the points is not in the list, it means it is farther than the farthest  point
+	 for (int j =0; j < NGB; j++)
+          {
+             if (&cells[i] == ngb[j]) {
+               found = true;
+               break;
+             }
+          } 
+	 double dist_to_seed = 0;
+	 for (int l = 0; l < 3; l++)
+	   {
+	     double delta = periodic(xc[k][l]-cells[i].coord[l]);
+	     dist_to_seed += delta*delta; 
+	   }
+	 double delta_epsilon = fabs(farther_dist/dist_to_seed-1);
+	 if (!found) 
+          { 
+	    if (dist_to_seed <= farther_dist &&
+		delta_epsilon > 1000*sqrt(numeric_limits<double>::epsilon()))
+	      {
+		cout << "SHOULD BE IN (d=" << dist_to_seed << "): Failed iteration " << k << " particle " << i << " (farthest=" << sqrt(farther_dist) << ")" << endl;
+		cout << "delta_epsilon = " << delta_epsilon << endl;
+		abort();
+	      }
+          }
+	 else
+	  {	    
+	    if (dist_to_seed > farther_dist && 
+		delta_epsilon > 1000*sqrt(numeric_limits<double>::epsilon())) 
+	      {	
+		cout << "SHOULD BE OUT (d=" << sqrt(dist_to_seed) << "): Failed iteration " << k << " particle " << i << " (farthest=" << sqrt(farther_dist) << ")" << endl;
+		cout << "delta_epsilon = " << delta_epsilon << " epsilon = " << 100*sqrt(numeric_limits<double>::epsilon()) << endl;
+		abort();
+	      }
+	  }
+      }
+  }
+
+  return 0;
+}

src/CMakeLists.txt

@@ -8,6 +8,7 @@ SET(CosmoTool_SRCS
 	powerSpectrum.cpp
 	miniargs.cpp
 	growthFactor.cpp
+	cosmopower.cpp
 )
 
 IF(FOUND_NETCDF3)

src/algo.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/algo.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMOTOOL_ALGO_HPP
 #define __COSMOTOOL_ALGO_HPP
 

src/bqueue.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/bqueue.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMO_QUEUE_HPP
 #define __COSMO_QUEUE_HPP
 

src/bsp_simple.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/bsp_simple.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMOTOOL_SIMPLE_BSP_HPP
 #define __COSMOTOOL_SIMPLE_BSP_HPP
 

src/cic.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/cic.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <assert.h>
 #include <math.h>
 #include <inttypes.h>

src/cic.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/cic.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __CICFILTER_HPP
 #define __CICFILTER_HPP
 

src/config.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/config.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMOTOOL_CONFIG_HPP
 #define __COSMOTOOL_CONFIG_HPP
 

src/cosmopower.cpp

@@ -0,0 +1,341 @@
+/*+
+This is CosmoTool (./src/cosmopower.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#include <cassert>
+#include <vector>
+#include <algorithm>
+#include <iostream>
+#include <cmath>
+#include <fstream>
+#include <gsl/gsl_integration.h>
+#include "cosmopower.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+#define USE_GSL
+#define TOLERANCE 1e-6
+#define NUM_ITERATION 8000
+
+
+CosmoPower::CosmoPower()
+{
+  eval = &CosmoPower::powerEfstathiou;
+  
+  n = 1.0;
+  K0 = 1;
+  V_LG_CMB = 627;
+
+  CMB_VECTOR[0] = 56.759;
+  CMB_VECTOR[1] = -540.02;
+  CMB_VECTOR[2] = 313.50;
+
+  h = 0.719;
+  SIGMA8 = 0.77;
+  OMEGA_B = 0.043969;
+  OMEGA_C = 0.21259;
+
+  Theta_27 = 2.728/2.7;
+
+  updateCosmology();
+}
+
+/*
+ * This is \hat{tophat}
+ */
+static double tophatFilter(double u)
+{
+  if (u != 0)    
+    return 3 / (u*u*u) * (sin(u) - u * cos(u));
+  else
+    return 1;
+}
+
+static double powC(double q, double alpha_c)
+{
+  return 14.2 / alpha_c + 386 / (1 + 69.9 * pow(q, 1.08));
+}
+
+static double T_tilde_0(double q, double alpha_c, double beta_c)
+{
+  double a = log(M_E + 1.8 * beta_c * q);
+  return a / ( a  + powC(q, alpha_c) * q * q);
+}
+
+static double j_0(double x)
+{
+  if (x == 0)
+    return 1.0;
+  return sin(x)/x;
+}
+
+static double powG(double y)
+{
+  double a = sqrt(1 + y);
+
+  return y * (-6 * a + (2 + 3 * y) *log((a + 1)/(a - 1)));
+}
+ 
+double CosmoPower::powerEfstathiou(double k)
+{
+  double a = 6.4/Gamma0;
+  double b = 3/Gamma0;
+  double c = 1.7/Gamma0;
+  double nu = 1.13;
+  
+  double f = (a*k) + pow(b*k,1.5) + pow(c*k,2);
+
+  // EFSTATHIOU  ET AL. (1992)
+  return normPower * pow(k,n) * pow(1+pow(f,nu),(-2/nu));
+}
+
+double CosmoPower::powerHuWiggles(double k)
+{
+  // EISENSTEIN ET HU (1998)
+  // FULL POWER SPECTRUM WITH BARYONS AND WIGGLES
+
+  double k_silk = 1.6 * pow(OMEGA_B * h * h, 0.52) * pow(OmegaEff, 0.73) * (1 + pow(10.4 * OmegaEff, -0.95));
+  double z_eq = 2.50e4 * OmegaEff * pow(Theta_27, -4);
+  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double f = 1 / (1 + pow(k * s / 5.4, 4));
+  double k_eq = 7.46e-2 * OmegaEff * pow(Theta_27, -2);
+  double a1 = pow(46.9 * OmegaEff, 0.670) * (1 + pow(32.1 * OmegaEff, -0.532));
+  double a2 = pow(12.0 * OmegaEff, 0.424) * (1 + pow(45.0 * OmegaEff, -0.582));
+  double alpha_c = pow(a1, -OMEGA_B/ OMEGA_0) * pow(a2, -pow(OMEGA_B / OMEGA_0, 3));  
+
+  double q = k / (13.41 * k_eq);
+  double b1_betac = 0.944 * 1/(1 + pow(458 * OmegaEff, -0.708));
+  double b2_betac = pow(0.395 * OmegaEff, -0.0266);
+  double beta_c = 1/ ( 1 + b1_betac * (pow(OMEGA_C / OMEGA_0, b2_betac) - 1)   );
+  double T_c = f * T_tilde_0(q, 1, beta_c) + (1 - f) * T_tilde_0(q, alpha_c, beta_c);
+
+  double b1_zd = 0.313 * pow(OmegaEff, -0.419) * (1 + 0.607 * pow(OmegaEff, 0.674));
+  double b2_zd = 0.238 * pow(OmegaEff, 0.223);
+  double z_d = 1291 * pow(OmegaEff, 0.251) / (1 + 0.659 * pow(OmegaEff, 0.828)) * (1 + b1_zd * pow(OmegaEff, b2_zd));
+  double R_d = 31.5 * OMEGA_B * h * h * pow(Theta_27, -4) * 1e3 / z_d;
+
+  double alpha_b = 2.07 * k_eq * s * pow(1 + R_d, -0.75) * powG((1 + z_eq)/(1 + z_d));
+  double beta_b = 0.5 + OMEGA_B / OMEGA_0 + (3 - 2 * OMEGA_B / OMEGA_0) * sqrt(pow(17.2 * OmegaEff, 2) + 1);
+  double beta_node = 8.41 * pow(OmegaEff, 0.435);
+  double s_tilde = s * pow(1 + pow(beta_node / (k * s), 3), -1./3);
+
+  double T_b = (T_tilde_0(q, 1, 1) / (1 + pow(k * s / 5.2, 2)) + alpha_b / (1 + pow(beta_b / (k * s), 3)) * exp(-pow(k/k_silk, 1.4))) * j_0(k * s_tilde);  
+
+  double T_k = OMEGA_B/OMEGA_0 * T_b + OMEGA_C/OMEGA_0 * T_c;  
+
+  return normPower * pow(k,n) * T_k * T_k;
+}
+
+double CosmoPower::powerHuBaryons(double k)
+{
+  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2);
+  double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4)));
+  double q = k/(h*GammaEff) * pow(Theta_27, 2);
+  double L_0 = log(2 * M_E + 1.8 * q);
+  double C_0 = 14.2 + 731 / (1 + 62.5 * q);
+  double T0 = L_0 / (L_0 + C_0 * q * q);
+  
+  return normPower * pow(k,n) * T0 * T0;
+}
+
+double CosmoPower::powerOld(double k)
+{
+  static const double l = 1/(Omega * h*h);
+  static const double alpha = 1.7 * l, beta = 9.0 * pow(l, 1.5), gamma = l*l;
+  return normPower * pow(k,n) * pow(1 + alpha * k +  beta * pow(k,1.5) + gamma *k*k,-2);
+}
+
+double CosmoPower::powerSugiyama(double k)
+{
+  double q = k * Theta_27*Theta_27 / (OmegaEff * exp(-OMEGA_B - sqrt(h/0.5)*OMEGA_B/OMEGA_0));
+  double L0 = log(2*M_E + 1.8 * q);
+  double C0 = 14.2 + 731 / (1 + 62.5 * q);
+  double T_k = L0 / (L0 + C0 * q*q);
+
+  return normPower * pow(k,n) * T_k * T_k;
+}
+
+double CosmoPower::powerBardeen(double k)
+{
+  double q = k / (OmegaEff);
+
+  double poly = 1 + 3.89 * q + pow(16.1*q,2) + pow(5.46*q,3) + pow(6.71*q,4);
+  double T_k = log(1+2.34*q)/(2.34*q) * pow(poly,-0.25);
+
+  return normPower * pow(k,n) * T_k * T_k;
+}
+
+double CosmoPower::powerBDM(double k)
+{
+  k /= h*h;
+  double alpha1 = 190;
+  double Gmu = 4.6;
+  double alpha2 = 11.5;
+  double alpha3 = 11;
+  double alpha4 = 12.55;
+  double alpha5 = 0.0004;
+  return normPower*k*alpha1*alpha1*Gmu*Gmu/(1+(alpha2*k)+pow(alpha3*k,2)+pow(alpha4*k,3))*pow(1+pow(alpha5/k,2), -2);
+}
+
+double CosmoPower::powerTest(double k)
+{
+  return 1/(1+k*k);
+}
+
+/*
+ * This function computes the normalization of the power spectrum. It requests
+ * a sigma8 (density fluctuations within 8 Mpc/h)
+ */
+static double gslPowSpecNorm(double k, void *params)
+{
+  CosmoPower *c = (CosmoPower *)params;
+
+  return c->integrandNormalize(k);
+}
+
+double CosmoPower::integrandNormalize(double x)
+{
+  double k = (1-x)/x;
+  double f = tophatFilter(k*8.0/h);
+  return power(k)*k*k*f*f/(x*x);
+}
+
+void CosmoPower::normalize()
+{
+  double normVal = 0;
+  double abserr;
+  gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION);
+  gsl_function f;
+
+  f.function = gslPowSpecNorm;
+  f.params = this;
+
+  normPower = 1;
+
+  ofstream ff("PP_k.txt");
+  for (int i = 0; i < 100; i++)
+    {
+       double k = pow(10.0, 4.0*i/100.-2);
+       ff << k << " " << power(k) << endl;
+    }
+
+ // gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr);
+  gsl_integration_qag(&f, 0, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr);
+  gsl_integration_workspace_free(w);
+
+  normVal /= (2*M_PI*M_PI);
+
+  normPower =  SIGMA8*SIGMA8/normVal;
+}
+
+void CosmoPower::updateCosmology()
+{
+  OMEGA_0 = OMEGA_B+OMEGA_C;
+  Omega = OMEGA_0;
+  beta = pow(OMEGA_0, 5./9);
+  OmegaEff = OMEGA_0 * h * h;
+  Gamma0 = OMEGA_0 * h * h;
+  omega_B = OMEGA_B * h * h;
+  omega_C = OMEGA_C * h * h;
+}
+
+void CosmoPower::updatePhysicalCosmology()
+{
+  OMEGA_B = omega_B / (h*h);
+  OMEGA_C = omega_C / (h*h);
+  OMEGA_0 = Gamma0 / (h*h);
+  beta = pow(OMEGA_0, 5./9);
+}
+
+double CosmoPower::eval_theta_theta(double k)
+{
+  // Jennings (2012) fit
+  double P_deltadelta = power(k);
+
+  static const double alpha0 = -12480.5, alpha1 = 1.824, alpha2 = 2165.87, alpha3=1.796;
+  if (k > 0.3)
+    return 0;
+  double r =(alpha0*sqrt(P_deltadelta) + alpha1*P_deltadelta*P_deltadelta)/(alpha2 + alpha3*P_deltadelta);
+  assert(P_deltadelta > 0);
+  
+  if (r < 0)
+    return 0;
+  return r;
+}
+
+double CosmoPower::power(double k)
+{
+  return (this->*eval)(k);
+}
+
+
+void CosmoPower::setFunction(CosmoFunction f)
+{
+  switch (f)
+    {
+    case POWER_EFSTATHIOU:
+      eval = &CosmoPower::powerEfstathiou;
+      break;
+    case HU_WIGGLES:
+      eval = &CosmoPower::powerHuWiggles;
+      break;
+    case HU_BARYON:
+      eval = &CosmoPower::powerHuBaryons;
+      break;
+    case OLD_POWERSPECTRUM:
+      eval = &CosmoPower::powerOld;
+      break;
+    case POWER_BARDEEN:
+      eval = &CosmoPower::powerBardeen;
+      break;
+    case POWER_SUGIYAMA:
+      eval = &CosmoPower::powerSugiyama;
+      break;
+    case POWER_BDM:
+      eval = &CosmoPower::powerBDM;
+      break;
+    case POWER_TEST:
+      eval = &CosmoPower::powerTest;
+      break;
+    default:
+      abort();
+    }
+}
+
+void CosmoPower::setNormalization(double A_K)
+{
+  normPower = A_K/power(0.002);
+}

src/cosmopower.hpp

@@ -0,0 +1,109 @@
+/*+
+This is CosmoTool (./src/cosmopower.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#ifndef _COSMOPOWER_HPP
+#define _COSMOPOWER_HPP
+
+namespace CosmoTool {
+
+  class CosmoPower
+  {
+  public:
+    // PRIMARY VARIABLES
+    double n;
+    double K0;
+    double V_LG_CMB;
+
+    double CMB_VECTOR[3];
+
+    double h;
+    double SIGMA8;
+    double OMEGA_B;
+    double OMEGA_C;
+    double omega_B;
+    double omega_C;
+    double Theta_27;
+    
+    // DERIVED VARIABLES    
+    double OMEGA_0;
+    double Omega;
+    double beta;
+    double OmegaEff;
+    double Gamma0;
+    double normPower;
+
+    enum CosmoFunction
+      {
+	POWER_EFSTATHIOU,
+	HU_WIGGLES,
+	HU_BARYON,
+	OLD_POWERSPECTRUM,
+	POWER_BARDEEN,
+	POWER_SUGIYAMA,
+	POWER_BDM,
+	POWER_TEST
+      };
+
+    CosmoPower();
+
+    void setFunction(CosmoFunction f);
+
+    void updateCosmology();
+    void updatePhysicalCosmology();
+    void normalize();
+    void setNormalization(double A_K);
+
+ 
+    double eval_theta_theta(double k);
+    double power(double k);
+
+    double integrandNormalize(double k);
+  private:
+    double (CosmoPower::*eval)(double);
+
+    double powerEfstathiou(double k);
+    double powerHuWiggles(double k);
+    double powerHuBaryons(double k);
+    double powerOld(double k);
+    double powerBardeen(double k);
+    double powerSugiyama(double k);
+    double powerBDM(double k);
+    double powerTest(double k);
+
+  };
+
+};
+
+#endif

src/dinterpolate.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/dinterpolate.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMO_DINTERPOLATE_HPP
 #define __COSMO_DINTERPOLATE_HPP
 
@@ -29,6 +64,7 @@ namespace CosmoTool {
     uint32_t numSimplex;
     uint32_t *simplex_list;
     gsl_eigen_symmv_workspace *eigen_work;
+    bool *disable_simplex;
     
     /**
      * This construct the interpolator. The construction is time consuming so
@@ -53,6 +89,7 @@ namespace CosmoTool {
       this->simplex_list = simplex_list;
       this->numPoints = numPoints;
       this->numSimplex = numSimplex;
+      this->disable_simplex = new bool[numSimplex];
       
       buildPreweight();
       buildQuickAccess();
@@ -67,6 +104,7 @@ namespace CosmoTool {
       delete quickAccess;
       delete[] point_to_simplex_list_base;      
       delete[] all_preweight;
+      delete[] disable_simplex;
 
       gsl_eigen_symmv_free(eigen_work);
     }

src/dinterpolate.tcc

@@ -1,3 +1,6 @@
+#include <fstream>
+#include <iostream>
+#include <iomanip>
 #include <cstdlib>
 #include <cassert>
 #include <gsl/gsl_matrix.h>
@@ -21,6 +24,7 @@ namespace CosmoTool {
     for (uint32_t i = 0; i < (N+1)*numSimplex; i++)
       {
 	assert(simplex_list[i] < numPoints);
+	if (!disable_simplex[i/(N+1)])
 	  numSimplex_by_point[simplex_list[i]]++;
       }
 
@@ -37,7 +41,11 @@ namespace CosmoTool {
       {
 	for (int j = 0; j <= N; j++)
 	  {
-	    uint32_t p = simplex_list[(N+1)*i+j];
+	    uint32_t s = (N+1)*i+j;
+	    if (disable_simplex[i])
+	      continue;
+
+	    uint32_t p = simplex_list[s];
 	    assert(index_by_point[p] < point_to_simplex_size);
 	    point_to_simplex_list_base[index_by_point[p]] = i;
 	    ++index_by_point[p];
@@ -48,7 +56,9 @@ namespace CosmoTool {
     for (uint32_t i = 0; i < numPoints; i++)
       {
 	// check assertion
-	assert((i==0 && index_by_point[0]==numSimplex_by_point[0]) || ((index_by_point[i]-index_by_point[i-1]) == (numSimplex_by_point[i]+1)));
+	assert((i==0 && index_by_point[0]==numSimplex_by_point[0]) 
+	       || 
+	       ((index_by_point[i]-index_by_point[i-1]) == (numSimplex_by_point[i]+1)));
 	assert(index_by_point[i] < point_to_simplex_size);
 	point_to_simplex_list_base[index_by_point[i]] = -1;
       }
@@ -80,6 +90,7 @@ namespace CosmoTool {
     double preweight[N*N];
     double preweight_inverse[N*N];
     gsl_permutation *p = gsl_permutation_alloc(N);
+    uint32_t numDisabled = 0;
 
     all_preweight = new PType[N*N*numSimplex];
 
@@ -105,13 +116,38 @@ namespace CosmoTool {
 	gsl_linalg_LU_decomp(&M.matrix, p, &signum);
         double a = fabs(gsl_linalg_LU_det(&M.matrix, signum));
 	if (a < 1e-10)
+         {
+#ifdef DEBUG
+           for (int j = 0; j < N; j++)
+            {
+              PType xref = positions[pref][j];
+
+              for (int k = 0; k < N; k++)
+                {
+                  preweight[j*N + k] = positions[simplex_list[k+base+1]][j] - xref;
+                }
+            }
+           std::ofstream f("matrix.txt");
+           for (int j = 0; j < N*N; j++)
+             f << std::setprecision(12) << preweight[j] << std::endl;
 	   throw InvalidArgumentException("Invalid tesselation. One simplex is coplanar.");
+#else
+           gsl_matrix_set_zero(&iM.matrix);
+	   disable_simplex[i] = true;
+	   numDisabled++;
+#endif
+         }
+        else {
 	  gsl_linalg_LU_invert(&M.matrix, p, &iM.matrix);
+	  disable_simplex[i] = false;
+        }
  
 	for (int j = 0; j < N*N; j++)
 	  all_preweight[N*N*i + j] = preweight_inverse[j];
       }
 
+    std::cout << "Number of disabled simplices: " << numDisabled << std::endl;
+
     gsl_permutation_free(p);
   }
 
@@ -166,6 +202,9 @@ namespace CosmoTool {
   template<typename PType, typename IType, int N>
   bool DelaunayInterpolate<PType,IType,N>::checkPointInSimplex(const CoordType& pos, uint32_t simplex)
   {
+    if (disable_simplex[simplex])
+      return false;
+
     uint32_t *desc_simplex = &simplex_list[simplex*(N+1)];
     CoordType *p[N+1], v[N], hyper;
 

src/eskow.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/eskow.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __ESKOW_CHOLESKY_HPP
 #define __ESKOW_CHOLESKY_HPP
 

src/field.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/field.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMOTOOL_FIELD
 #define __COSMOTOOL_FIELD
 

src/fixArray.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/fixArray.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __FIX_ARRAY_HPP
 #define __FIX_ARRAY_HPP
 

src/fortran.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/fortran.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include "config.hpp"
 #include <iostream>
 #include <fstream>

src/fortran.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/fortran.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMO_FORTRAN_HPP
 #define __COSMO_FORTRAN_HPP
 

src/fourier/base_types.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/fourier/base_types.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __BASE_FOURIER_TYPES_HPP
 #define __BASE_FOURIER_TYPES_HPP
 
@@ -24,14 +59,28 @@ namespace CosmoTool
     typedef Eigen::Map<VecType, Eigen::Aligned> MapType;
     typedef Eigen::Map<VecType const, Eigen::Aligned> ConstMapType;
     typedef FourierMap<std::complex<T> > FourierMapType;
+    typedef boost::shared_ptr<FourierMapType> FourierMapPtr;
+    typedef boost::shared_ptr<SpectrumFunction<T> > SpectrumFunctionPtr;
 
-    virtual boost::shared_ptr<FourierMapType> 
+    virtual ~SpectrumFunction() {}
-       newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0;
 
-    virtual void mul(FourierMap<std::complex<T> >& m) const = 0;
+    virtual void 
+       newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const = 0;
+
+    virtual SpectrumFunctionPtr copy() const = 0;
+
+    virtual const T *data() const { return 0; }
+    virtual const int *dof() const { return 0; }
+    virtual long size() const { return -1; }
+
+    virtual void sqrt() = 0;
+
+    virtual void mul(FourierMapType& m) const = 0;
+    virtual void mul_sqrt(FourierMapType& m) const = 0;
+    virtual void mul_inv(FourierMapType& m) const = 0;
+    virtual void mul_inv_sqrt(FourierMapType& m) const = 0;
   };
 
-
   template<typename T>
   class FourierMap
   {
@@ -60,6 +109,12 @@ namespace CosmoTool
       return ConstMapType(data(), size());
     }
 
+    FourierMap<T>& operator=(const FourierMap<T>& m)
+    {
+      eigen() = m.eigen();
+      return *this;
+    }
+
     void sqrt()
     {
       MapType m = eigen();
@@ -76,7 +131,7 @@ namespace CosmoTool
     {
       assert(size() == map2->size());
       MapType m(data(), size());
-      MapType m2(map2->data(), map2->size());
+      ConstMapType m2(map2->data(), map2->size());
       m *= m2;
     }
 
@@ -130,6 +185,18 @@ namespace CosmoTool
     virtual void synthesis_conjugate() = 0;
   };
 
+  template<typename T>
+  class MapUtilityFunction
+  {
+  public:
+    typedef SpectrumFunction<T> Spectrum;
+    typedef boost::shared_ptr<Spectrum> Spectrum_ptr;
+    typedef FourierMap<std::complex<T> > FMap;
+
+    virtual Spectrum_ptr estimateSpectrumFromMap(const FMap& m) const = 0;
+    virtual Spectrum_ptr newSpectrumFromRaw(T *data, long size,
+                                            const Spectrum& like_spec) const = 0;
+  };
 
 };
 

src/fourier/details/euclidian_maps.hpp

@@ -0,0 +1,267 @@
+/*+
+This is CosmoTool (./src/fourier/details/euclidian_maps.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#ifndef __DETAILS_EUCLIDIAN_MAPS
+#define __DETAILS_EUCLIDIAN_MAPS
+
+
+namespace CosmoTool
+{
+  
+  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;
+    int m_dim0;
+    bool even0, alleven;
+    long plane_size;
+  public:
+    typedef typename EuclidianFourierMapBase<std::complex<T> >::DimArray DimArray;
+
+    EuclidianFourierMapComplex(ptr_t indata,
+			       int dim0,
+			       const DimArray& indims, 
+			       const std::vector<double>& dk)
+      : EuclidianFourierMapBase<std::complex<T> >(indata, indims), delta_k(dk), m_dim0(dim0), even0((dim0 % 2)==0)
+    {
+      assert(dk.size() == indims.size()); 
+      plane_size = 1;
+      alleven = true;
+      for (int q = 1; q < indims.size(); q++)
+	{
+	  plane_size *= indims[q];
+	  alleven = alleven && ((indims[q]%2)==0);
+	}
+    }
+
+    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)),
+		   m_dim0,
+		   this->getDims(),
+		   this->delta_k);
+    }
+
+    const std::vector<double>& get_delta_k() const
+    {
+      return this->delta_k;
+    }
+
+    template<typename Array, typename Array2>
+    void get_Kvec(const Array& ik, Array2& kvec) const
+    {
+      const DimArray& dims = this->getDims();
+      assert(ik.size() == dims.size());
+      assert(kvec.size() == dims.size());
+
+      kvec[0] = ik[0] * delta_k[0];
+      for (int q = 1; q < ik.size(); q++)
+	{
+	  int dk = ik[q];
+	  if (dk > dims[q]/2)
+	    dk = dk - dims[q];
+
+          kvec[q] = dk*delta_k[q];
+        }
+    }
+
+    template<typename Array2>
+    void get_Kvec(long p, Array2& kvec) const
+    {
+      const DimArray& dims = this->getDims();
+      DimArray d(delta_k.size());
+      get_IKvec(p, d);
+      get_Kvec(d, kvec);
+    }
+
+    void get_IKvec(long p, DimArray& ikvec) const
+    {
+      const DimArray& dims = this->getDims();
+      assert(dims.size()==ikvec.size());
+      for (int q = 0; q < ikvec.size(); q++)
+       {
+         ikvec[q] = p%dims[q];
+         p = (p-ikvec[q])/dims[q];
+       }
+    }
+
+
+    template<typename Array>
+    double get_K(const Array& ik) const
+    {
+      const DimArray& dims = this->getDims();
+      assert(ik.size() == dims.size());
+      double k2 = 0;
+      k2 += CosmoTool::square(ik[0]*delta_k[0]);
+
+      for (int q = 1; q < ik.size(); q++)
+	{
+	  int dk = ik[q];
+
+	  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
+    {
+      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);
+    }
+
+    bool allDimensionsEven() const { return alleven; }
+    bool firstDimensionEven() const { return even0; }
+
+    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] + (even0 ? 0 : 1);
+      std::complex<T> result = 0;
+
+      for (long q0 = 1; q0 < N0-1; q0++)
+	{
+	  for (long p = 0; p < plane_size; p++)
+	    {
+	      long idx = q0+dims[0]*p;
+	      assert(idx < this->size());
+	      result += 2*(conj(d1[idx]) * d2[idx]).real();
+	    }
+	}
+      if (even0)
+	{
+	  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;
+    }
+
+  };
+};
+
+#endif

src/fourier/details/euclidian_spectrum_1d.hpp

@@ -0,0 +1,215 @@
+/*+
+This is CosmoTool (./src/fourier/details/euclidian_spectrum_1d.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#ifndef __DETAILS_EUCLIDIAN_SPECTRUM_1D
+#define __DETAILS_EUCLIDIAN_SPECTRUM_1D
+
+#include  <iostream>
+#include <boost/function.hpp>
+
+
+namespace CosmoTool
+{
+  template<typename T>
+  class EuclidianOperator
+  {
+  public:
+    typedef boost::function1<T, T> Function;
+    
+    Function base, op;
+    T operator()(T k) {
+      return op(base(k));
+    }
+  };
+
+  template<typename T>
+  class EuclidianSpectrum_1D: public SpectrumFunction<T>
+  {
+  public:
+    typedef boost::function1<T, T> Function;
+  protected:
+    Function f;
+
+    static T msqrt(T a) { return std::sqrt(a); }
+  public:
+    typedef typename SpectrumFunction<T>::FourierMapType FourierMapType;
+    typedef typename SpectrumFunction<T>::SpectrumFunctionPtr SpectrumFunctionPtr;
+    typedef boost::shared_ptr<FourierMapType> ptr_map;
+
+    EuclidianSpectrum_1D(Function P)
+      : f(P)
+    {
+    }
+
+    void newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const;
+
+    SpectrumFunctionPtr copy() const {
+      return SpectrumFunctionPtr(new EuclidianSpectrum_1D(f));
+    }
+
+    void sqrt() {
+      EuclidianOperator<T> o;
+      o.base = f;
+      o.op = &EuclidianSpectrum_1D<T>::msqrt;
+      f = (Function(o));
+    }
+
+    void mul(FourierMapType& m) const;
+    void mul_sqrt(FourierMapType& m) const;
+    void mul_inv(FourierMapType& m) const;
+    void mul_inv_sqrt(FourierMapType& m) const;
+  };
+
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const
+  {
+    typedef EuclidianFourierMapComplex<T> MapT;
+    typedef typename EuclidianSpectrum_1D<T>::ptr_map ptr_map;
+    typedef typename MapT::DimArray DimArray;
+
+    MapT& rand_map = dynamic_cast<MapT&>(out_map);
+
+    std::complex<T> *d = rand_map.data();
+    long idx;
+    const DimArray& dims = rand_map.getDims();
+    const std::vector<double>& delta_k = rand_map.get_delta_k();
+    long plane_size;
+    bool alleven = rand_map.allDimensionsEven();
+    double V = 1;
+
+    for (int p = 0; p < delta_k.size(); p++)
+      V *= (2*M_PI/delta_k[p]);
+
+    for (long p = 1; p < rand_map.size(); p++)
+      {
+	double A_k = std::sqrt(0.5*V*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>(gsl_ran_gaussian(rng, std::sqrt(V*f(0))), 0);
+
+    if (!rand_map.firstDimensionEven())
+      return;
+
+    // Correct the Nyquist plane
+    idx = dims[0]-1; // Stick to the last element of the first dimension
+    d[idx] = std::complex<T>(d[idx].real() + d[idx].imag(), 0);
+    // 1D is special case
+    if (dims.size() == 1)
+      return;
+    
+    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]);
+      }
+
+    if (alleven)
+      {
+	long q = 0;
+	for (int i = dims.size()-1; i >= 1; i--)
+	  q = dims[i]*q + dims[i]/2;
+	q += dims[0]-1;
+	d[q] = std::complex<T>(d[q].real()+d[q].imag(),0);
+      }
+  }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul(FourierMapType& 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_c.get_K(p));
+  }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul_sqrt(FourierMapType& 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] *= std::sqrt(f(m_c.get_K(p)));
+  }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul_inv(FourierMapType& 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++)
+     {
+        T A = f(m_c.get_K(p));
+        if (A==0)
+          d[p] = 0;
+        else
+          d[p] /= A;
+     }
+  }
+
+  template<typename T>
+  void EuclidianSpectrum_1D<T>::mul_inv_sqrt(FourierMapType& 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++)
+      {
+        T A = std::sqrt(f(m_c.get_K(p)));
+        if (A == 0)
+          d[p] = 0;
+        else
+          d[p] /= A;
+      }
+  }
+
+};
+
+
+#endif

src/fourier/details/euclidian_spectrum_1d_bin.hpp

@@ -0,0 +1,100 @@
+/*+
+This is CosmoTool (./src/fourier/details/euclidian_spectrum_1d_bin.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#ifndef __DETAILS_EUCLIDIAN_SPECTRUM_1D_BIN
+#define __DETAILS_EUCLIDIAN_SPECTRUM_1D_BIN
+
+#include <boost/bind.hpp>
+#include <cmath>
+
+namespace CosmoTool
+{
+  template<typename T>
+  class EuclidianSpectrum_1D_Binned: public EuclidianSpectrum_1D<T>
+  {
+  protected:
+    T *m_data;
+    long m_size;
+    T m_kmin, m_kmax, m_logdeltak;
+    std::vector<T> m_dof;
+  public:
+    typedef typename SpectrumFunction<T>::FourierMapType FourierMapType;
+    typedef typename SpectrumFunction<T>::SpectrumFunctionPtr SpectrumFunctionPtr;
+    typedef boost::shared_ptr<FourierMapType> ptr_map;
+
+    T interpolate_spectrum(T k)
+    {
+      T q = std::log(k/m_kmin)/m_logdeltak;
+      long ik = std::floor(q);
+
+      if (ik >= m_size-1)
+	return m_data[m_size-1];
+      else if (ik < 0)
+	return k/m_kmin*m_data[0];
+
+      return std::exp((q-ik)*m_data[ik+1] + (1-(q-ik))*m_data[ik]);
+    }
+
+    EuclidianSpectrum_1D_Binned(int numBin, T kmin, T kmax)
+      : EuclidianSpectrum_1D<T>(boost::bind(&EuclidianSpectrum_1D_Binned::interpolate_spectrum, this, _1))
+    {
+      m_data = new T[numBin];
+      m_kmin = kmin;
+      m_kmax = kmax;
+      m_size = numBin;
+      m_logdeltak = std::log(m_kmax/m_kmin);
+    }
+
+    SpectrumFunctionPtr copy() const
+    {
+      EuclidianSpectrum_1D_Binned *s = new EuclidianSpectrum_1D_Binned(m_size, m_kmin, m_kmax);
+      std::copy(m_data, m_data+m_size, s->m_data);
+      return SpectrumFunctionPtr(s);
+    }
+
+    void set_dof(std::vector<T>& dof_array) 
+    {
+      assert(m_size == dof_array.size());
+      m_dof = dof_array;
+    }
+
+    const T *data() const { return m_data; }
+    long size() const { return m_size; }
+    const T *dof() const { return &m_dof[0]; }
+  };
+
+};
+
+#endif

src/fourier/details/euclidian_transform.hpp

@@ -0,0 +1,216 @@
+/*+
+This is CosmoTool (./src/fourier/details/euclidian_transform.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#ifndef __DETAILS_EUCLIDIAN_TRANSFORM
+#define __DETAILS_EUCLIDIAN_TRANSFORM
+
+namespace CosmoTool
+{
+  
+  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, Nc;
+    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());
+      std::vector<int> swapped_dims(dims.size());
+      
+      m_dims = dims;
+      m_dims_hc = dims;
+      m_dims_hc[0] = dims[0]/2+1;
+      m_L = L;
+
+      N = 1;
+      Nc = 1;
+      volume = 1;
+      for (int i = 0; i < dims.size(); i++)
+	{
+	  N *= dims[i];
+          Nc *= m_dims_hc[i];
+	  volume *= L[i];
+	  dk[i] = 2*M_PI/L[i];
+          swapped_dims[dims.size()-1-i] = dims[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(Nc),
+						       std::ptr_fun(calls::free)), 
+		   dims[0], m_dims_hc, dk);
+      m_analysis = calls::plan_dft_r2c(dims.size(), &swapped_dims[0],
+				       realMap->data(), (typename calls::complex_type *)fourierMap->data(),
+				       FFTW_DESTROY_INPUT|FFTW_MEASURE);
+      m_synthesis = calls::plan_dft_c2r(dims.size(), &swapped_dims[0],
+					(typename calls::complex_type *)fourierMap->data(), realMap->data(),
+					FFTW_DESTROY_INPUT|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_1d: public EuclidianFourierTransform<T>
+  {
+  private:
+    template<typename T2>
+    static std::vector<T2> make_1d_vector(T2 a)
+    {
+      T2 arr[2] = { a};
+      return std::vector<T2>(&arr[0],&arr[1]);
+    }
+  public:
+    EuclidianFourierTransform_1d(int Nx, double Lx)
+      : EuclidianFourierTransform<T>(make_1d_vector<int>(Nx), make_1d_vector<double>(Lx))
+    {
+    }
+
+    virtual ~EuclidianFourierTransform_1d() {}
+
+  };
+
+  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[3] = { a, b, c};
+      return std::vector<T2>(&arr[0], &arr[3]);
+    }
+
+  public:
+    EuclidianFourierTransform_3d(int Nx, int Ny, int Nz, double Lx, double Ly, double Lz)
+      : EuclidianFourierTransform<T>(make_3d_vector<int>(Nx, Ny, Nz), make_3d_vector<double>(Lx, Ly, Lz))
+    {
+    }
+
+    virtual ~EuclidianFourierTransform_3d() {}
+  };
+
+};
+
+#endif

src/fourier/details/healpix_alms.hpp

@@ -0,0 +1,111 @@
+/*+
+This is CosmoTool (./src/fourier/details/healpix_alms.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_ALM_HPP
+#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_ALM_HPP
+
+namespace CosmoTool
+{
+  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;
+    }
+  };
+};
+
+#endif

src/fourier/details/healpix_map.hpp

@@ -0,0 +1,89 @@
+/*+
+This is CosmoTool (./src/fourier/details/healpix_map.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_MAP_HPP
+#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_MAP_HPP
+
+namespace CosmoTool
+{
+
+  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);
+    }
+  };
+
+};
+
+#endif

src/fourier/details/healpix_spectrum.hpp

@@ -0,0 +1,179 @@
+/*+
+This is CosmoTool (./src/fourier/details/healpix_spectrum.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HPP
+#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HPP
+
+namespace CosmoTool
+{
+  template<typename T>
+  class HealpixSpectrum: public SpectrumFunction<T>
+  {
+  protected:
+    std::vector<T> cls;
+    int *m_dof;
+  public:
+    typedef typename SpectrumFunction<T>::FourierMapType FourierMapType;
+    typedef boost::shared_ptr<FourierMapType> ptr_map;
+    typedef typename SpectrumFunction<T>::SpectrumFunctionPtr SpectrumFunctionPtr;
+    typedef unsigned long LType;
+
+    HealpixSpectrum(LType Lmax)
+      : cls(Lmax+1), m_dof(new int[Lmax+1]) 
+    {
+      for (int l = 0; l <= Lmax; l++)
+	m_dof[l] = l + 1;
+    }
+
+    T *data() { return &cls[0]; }
+    const T *data() const { return &cls[0]; }
+    const int *dof() const { return m_dof; } 
+    
+    void set_dof(LType l, int dof) { m_dof[l] = dof; }
+
+    LType Lmax() const { return cls.size()-1; }
+    long size() const { return cls.size(); }
+
+    void newRandomFourier(gsl_rng *rng, FourierMapType& like_map) const; 
+
+    SpectrumFunctionPtr copy() const {
+      HealpixSpectrum<T> *s = new HealpixSpectrum<T>(Lmax());
+      s->cls = cls;
+      return SpectrumFunctionPtr(s);
+    }
+
+    void sqrt() {
+      std::transform(cls.begin(), cls.end(), cls.begin(), std::ptr_fun<T,T>(std::sqrt));
+    }
+
+    void mul(FourierMapType& m) const;
+    void mul_sqrt(FourierMapType& m) const;
+    void mul_inv(FourierMapType& m) const;
+    void mul_inv_sqrt(FourierMapType& m) const;
+  };
+
+  template<typename T>
+  void HealpixSpectrum<T>::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const
+  { 
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(out_map);
+    long lmaxGen = std::min(cls.size()-1, alms.Lmax());
+    std::complex<T> *new_data = alms.data();
+
+    for (LType 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 (LType m = 1; m <= std::min(l,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); 
+          }
+      }
+  } 
+
+ template<typename T>
+ void HealpixSpectrum<T>::mul(FourierMapType& like_map) const
+  {
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
+    std::complex<T> *data = alms.data();
+
+    for (LType l = 0; l <= alms.Lmax(); l++)
+      {
+        double Al = cls[l];
+
+        for (LType m = 0; m <= std::min(l,alms.Mmax()); m++)
+          {
+            data[alms.index(l,m)] *= Al;
+          }
+      }
+  }
+
+ template<typename T>
+ void HealpixSpectrum<T>::mul_sqrt(FourierMapType& like_map) const
+  {
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
+    std::complex<T> *data = alms.data();
+
+    for (LType l = 0; l <= alms.Lmax(); l++)
+      {
+        double Al = std::sqrt(cls[l]);
+
+        for (LType m = 0; m <= std::min(l,alms.Mmax()); m++)
+          {
+            data[alms.index(l,m)] *= Al;
+          }
+      }
+  }
+
+ template<typename T>
+ void HealpixSpectrum<T>::mul_inv(FourierMapType& like_map) const
+  {
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
+    std::complex<T> *data = alms.data();
+
+    for (LType l = 0; l <= alms.Lmax(); l++)
+      {
+        double Al = (cls[l] <= 0) ? 0 : (1/cls[l]);
+
+        for (LType m = 0; m <= std::min(l,alms.Mmax()); m++)
+          {
+            data[alms.index(l,m)] *= Al;
+          }
+      }
+  }
+
+ template<typename T>
+ void HealpixSpectrum<T>::mul_inv_sqrt(FourierMapType& like_map) const
+  {
+    HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(like_map);
+    std::complex<T> *data = alms.data();
+
+    for (LType l = 0; l <= alms.Lmax(); l++)
+      {
+        double Al = (cls[l] <= 0) ? 0 : std::sqrt(1/cls[l]);
+
+        for (LType m = 0; m <= std::min(l,alms.Mmax()); m++)
+          {
+            data[alms.index(l,m)] *= Al;
+          }
+      }
+  }
+
+};
+
+#endif

src/fourier/details/healpix_transform.hpp

@@ -0,0 +1,132 @@
+/*+
+This is CosmoTool (./src/fourier/details/healpix_transform.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_TRANSFORM_HPP
+#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_TRANSFORM_HPP
+
+namespace CosmoTool
+{
+
+  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));
+    }
+    
+  };
+
+};
+
+#endif

src/fourier/details/healpix_utility.hpp

@@ -0,0 +1,98 @@
+/*+
+This is CosmoTool (./src/fourier/details/healpix_utility.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HP
+#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HP
+
+namespace CosmoTool
+{
+
+  template<typename T>
+  class HealpixUtilityFunction: public MapUtilityFunction<T>
+  {
+  public:
+    typedef typename MapUtilityFunction<T>::Spectrum Spectrum;
+    typedef typename MapUtilityFunction<T>::Spectrum_ptr Spectrum_ptr;
+    typedef typename MapUtilityFunction<T>::FMap FMap;
+    typedef typename HealpixSpectrum<T>::LType LType; 
+    
+    Spectrum_ptr estimateSpectrumFromMap(const FMap& m) const
+    {
+      const HealpixFourierALM<T>& alms = dynamic_cast<const HealpixFourierALM<T>&>(m);
+      LType Lmax = alms.Lmax(), Mmax = alms.Mmax();
+      HealpixSpectrum<T> *spectrum = new HealpixSpectrum<T>(Lmax);
+      T *data_cls = spectrum->data();
+      const std::complex<T> *data_alms = alms.data();
+
+      // make an estimate of the cls
+      std::fill(data_cls, data_cls+Lmax+1, 0);
+      LType q = 0;
+      for (LType m = 0; m <= Mmax; ++m )
+	{
+	  for (LType l = m; l <= Lmax; ++l)
+	    {
+	      // Triangular storage
+	      data_cls[l] += std::norm(data_alms[l+q]);
+	    }
+	  q += Lmax+1-m;
+	}
+      assert(q==alms.size());
+      
+      for (LType l = 0; l <= Lmax; ++l)
+	{
+	  int dof = 1 + std::min(l, Mmax);
+	  spectrum->set_dof(l, dof);
+	  data_cls[l] /= dof;
+	}
+
+      return Spectrum_ptr(spectrum);
+    }
+
+    Spectrum_ptr newSpectrumFromRaw(T *data, long size,
+				    const Spectrum& like_spec) const
+    {
+      const HealpixSpectrum<T>& in_spec = dynamic_cast<const HealpixSpectrum<T>&>(like_spec);
+      HealpixSpectrum<T> *new_spectrum = new HealpixSpectrum<T>(in_spec.Lmax());
+      T *out_d = new_spectrum->data();
+
+      std::copy(data, data + min(size,new_spectrum->size()), out_d);
+
+      return Spectrum_ptr(new_spectrum);
+    }
+  };
+
+};
+
+#endif

src/fourier/euclidian.hpp

@@ -1,6 +1,42 @@
+/*+
+This is CosmoTool (./src/fourier/euclidian.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
 #define __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
 
+#include <cmath>
 #include <boost/function.hpp>
 #include <vector>
 #include <boost/shared_ptr.hpp>
@@ -8,400 +44,9 @@
 #include "base_types.hpp"
 #include "fft/fftw_calls.hpp"
 #include "../algo.hpp"
-
+#include "details/euclidian_maps.hpp"
-namespace CosmoTool
+#include "details/euclidian_spectrum_1d.hpp"
-{
+#include "details/euclidian_spectrum_1d_bin.hpp"
-  template<typename T>
+#include "details/euclidian_transform.hpp"
-  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

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/fourier/fft/fftw_calls.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __FFTW_UNIFIED_CALLS_HPP
 #define __FFTW_UNIFIED_CALLS_HPP
 
@@ -28,8 +63,8 @@ public: \
   typedef prefix ## _complex complex_type; \
   typedef prefix ## _plan plan_type; \
   \
-  static complex_type *alloc_complex(int N) { return prefix ## _alloc_complex(N); } \
+  static complex_type *alloc_complex(size_t N) { return prefix ## _alloc_complex(N); } \
-  static real_type *alloc_real(int N) { return prefix ## _alloc_real(N); } \
+  static real_type *alloc_real(size_t N) { return prefix ## _alloc_real(N); } \
   static void free(void *p) { fftw_free(p); } \
 \
   static void execute(plan_type p) { prefix ## _execute(p); } \

src/fourier/healpix.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/fourier/healpix.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMOTOOL_FOURIER_HEALPIX_HPP
 #define __COSMOTOOL_FOURIER_HEALPIX_HPP
 
@@ -12,257 +47,13 @@
 #include <sharp_lowlevel.h>
 #include <sharp_geomhelpers.h>
 #include <sharp_almhelpers.h>
+#include <algorithm>
+#include <functional>
 
-namespace CosmoTool
+#include "details/healpix_map.hpp"
-{
+#include "details/healpix_alms.hpp"
-  template<typename T>
+#include "details/healpix_transform.hpp"
-  class HealpixSpectrum: public SpectrumFunction<T>
+#include "details/healpix_spectrum.hpp"
-  {
+#include "details/healpix_utility.hpp"
-  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/growthFactor.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/growthFactor.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <cmath>
 #include <gsl/gsl_integration.h>
 #include "interpolate.hpp"

src/growthFactor.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/growthFactor.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef COSMO_GROWTH_FACTOR_HPP
 #define COSMO_GROWTH_FACTOR_HPP
 

src/h5_readFlash.cpp

@@ -1,3 +1,11 @@
+/*+
+!!
+
+This file has been developped by P. M. Sutter
+
+!!
+
++*/
 /* This file contains the functions that read the data from the HDF5 file
  * The functions accept the PARAMESH data through arguments.
  */
@@ -185,7 +193,7 @@ void h5_read_flash3_particles (H5File* file,
                         float *vel1,
                         float *vel2,
                         float *vel3,
-                        int    *id)
+                        long    *id)
 {
 
   herr_t   status;
@@ -341,7 +349,7 @@ void h5_read_flash3_particles (H5File* file,
 
     if (id) {
     for(p=0; p < (pcount); p++) {
-      id[p+poffset] = (int)  *(partBuffer+iptag-1+p*numProps);
+      id[p+poffset] = (long)  *(partBuffer+iptag-1+p*numProps);
     } }
 
     if (pos1 && pos2 && pos3) {

src/h5_readFlash.hpp

@@ -1,3 +1,12 @@
+/*+
+
+!!! NOTE !!!
+
+This file has been developped by P. M. Sutter.
+
+!!!!
+
++*/
 /* This file contains the functions that read the data from the HDF5 file
  * The functions accept the PARAMESH data through arguments.
  */
@@ -29,7 +38,7 @@ void h5_read_flash3_particles (H5File* file,
                         float *vel1,
                         float *vel2,
                         float *vel3,
-                        int    *id);
+                        long    *id);
 
 void h5_read_flash3_header_info(H5File* file,
 				double* time,                   /* simulation time */

src/hdf5_flash.h

@@ -1,3 +1,10 @@
+/*+
+!!
+
+This particular file has been developped by P. M. Sutter
+
+!!
++*/
 /* general header file for the HDF 5 IO in FLASH */
 
 

src/interpolate.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/interpolate.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <cassert>
 #include <cmath>
 #include <algorithm>

src/interpolate.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/interpolate.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __CTOOL_INTERPOLATE_HPP
 #define __CTOOL_INTERPOLATE_HPP
 

src/interpolate3d.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/interpolate3d.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMO_INTERPOLATE3D_HPP
 #define __COSMO_INTERPOLATE3D_HPP
 

src/kdtree_leaf.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/kdtree_leaf.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __LEAF_KDTREE_HPP
 #define __LEAF_KDTREE_HPP
 

src/kdtree_splitters.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/kdtree_splitters.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __KDTREE_SPLITTERS_HPP
 #define __KDTREE_SPLITTERS_HPP
 

src/loadFlash.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/loadFlash.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 /* Reads in FLASH v3 files in HDF5 format */
 
 #include <iostream>
@@ -71,7 +106,7 @@ SimuData *CosmoTool::loadFlashMulti(const char *fname, int id, int loadflags)
     } }
  
     if (loadflags & NEED_GADGET_ID) {
-    data->Id = new int[data->NumPart];
+    data->Id = new long[data->NumPart];
     if (data->Id == 0) {
       delete data;
       return 0;

src/loadFlash.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/loadFlash.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMO_LOAD_FLASH_HPP
 #define __COSMO_LOAD_FLASH_HPP
 

src/loadFlash_dummy.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/loadFlash_dummy.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include "load_data.hpp"
 #include "loadFlash.hpp"
 

src/loadGadget.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/loadGadget.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <cmath>
 #include <iostream>
 #include <assert.h>
@@ -10,61 +45,41 @@
 using namespace CosmoTool;
 using namespace std;
 
-PurePositionData *CosmoTool::loadGadgetPosition(const char *fname)
+
+void loadGadgetHeader(UnformattedRead *f, GadgetHeader& h, SimuData *data, int id)
 {
-  PurePositionData *data;
+  f->beginCheckpoint();
-  int p, n;
-  UnformattedRead f(fname);
-  GadgetHeader h;
-  
-  data = new PurePositionData;
-  f.beginCheckpoint();
   for (int i = 0; i < 6; i++)
-    h.npart[i] = f.readInt32();
+    h.npart[i] = f->readInt32();
   for (int i = 0; i < 6; i++)    
-    h.mass[i] = f.readReal64();
+    h.mass[i] = f->readReal64();
-  h.time = f.readReal64();
+  data->time = h.time = f->readReal64();
-  h.redshift = f.readReal64();
+  h.redshift = f->readReal64();
-  h.flag_sfr = f.readInt32();
+  h.flag_sfr = f->readInt32();
-  h.flag_feedback = f.readInt32();
+  h.flag_feedback = f->readInt32();
   for (int i = 0; i < 6; i++)
-    h.npartTotal[i] = f.readInt32();
+    h.npartTotal[i] = f->readInt32();
-  h.flag_cooling = f.readInt32();
+  h.flag_cooling = f->readInt32();
-  h.num_files = f.readInt32();
+  h.num_files = f->readInt32();
-  data->BoxSize = h.BoxSize = f.readReal64();
+  data->BoxSize = h.BoxSize = f->readReal64();
-  h.Omega0 = f.readReal64();
+  data->Omega_M = h.Omega0 = f->readReal64();
-  h.OmegaLambda = f.readReal64();
+  data->Omega_Lambda = h.OmegaLambda = f->readReal64();
-  h.HubbleParam = f.readReal64();
+  data->Hubble = h.HubbleParam = f->readReal64();
-  f.endCheckpoint(true);
+  f->endCheckpoint(true);
   
-  data->NumPart = 0;
+  long NumPart = 0, NumPartTotal = 0;
-  for(int k=0; k<5; k++)
+  for(int k=0; k<6; k++)
-    data->NumPart += h.npart[k];
+    {
-  
+      NumPart += h.npart[k];
-  data->pos = new FCoordinates[data->NumPart];
+      NumPartTotal += (id < 0) ? h.npart[k] : h.npartTotal[k];
-  
-  f.beginCheckpoint();
-  for(int k = 0, p = 0; k < 5; k++) {
-    for(int n = 0; n < h.npart[k]; n++) {
-      data->pos[p][0] = f.readReal32();
-      data->pos[p][1] = f.readReal32();
-      data->pos[p][2] = f.readReal32();
-      p++;
     }
-  }
+  data->NumPart = NumPart;
-  f.endCheckpoint();
+  data->TotalNumPart = NumPartTotal;
-  
-  // Skip velocities
-  f.skip((long)data->NumPart*3+2*4);
-  // Skip ids
-
-  return data;
 }
 
-
+SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
-
+				     int loadflags, int GadgetFormat, 
-
+				     SimuFilter filter)
-SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, int GadgetFormat)
 {
   SimuData *data;
   int p, n;
@@ -101,32 +116,8 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i
   
   try
     {
-      f->beginCheckpoint();
+      loadGadgetHeader(f, h, data, id);
-      for (int i = 0; i < 6; i++)
-	h.npart[i] = f->readInt32();
-      for (int i = 0; i < 6; i++)    
-	h.mass[i] = f->readReal64();
-      data->time = h.time = f->readReal64();
-      h.redshift = f->readReal64();
-      h.flag_sfr = f->readInt32();
-      h.flag_feedback = f->readInt32();
-      for (int i = 0; i < 6; i++)
-	h.npartTotal[i] = f->readInt32();
-      h.flag_cooling = f->readInt32();
-      h.num_files = f->readInt32();
-      data->BoxSize = h.BoxSize = f->readReal64();
-      data->Omega_M = h.Omega0 = f->readReal64();
-      data->Omega_Lambda = h.OmegaLambda = f->readReal64();
-      data->Hubble = h.HubbleParam = f->readReal64();
-      f->endCheckpoint(true);
 
-      for(int k=0; k<6; k++)
-	{
-	  NumPart += h.npart[k];
-	  NumPartTotal += (id < 0) ? h.npart[k] : h.npartTotal[k];
-	}
-      data->NumPart = NumPart;
-      data->TotalNumPart = NumPartTotal;
       if (GadgetFormat == 1)
         velmul = sqrt(h.time);
       else if (GadgetFormat == 2)
@@ -135,6 +126,9 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i
         cerr << "unknown gadget format" << endl;
 	abort();
       }
+
+      NumPart = data->NumPart;
+      NumPartTotal = data->TotalNumPart;
     }
   catch (const InvalidUnformattedAccess& e)
     {
@@ -236,7 +230,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i
     try
       {
 	f->beginCheckpoint();
-	data->Id = new int[data->NumPart];
+	data->Id = new long[data->NumPart];
 	if (data->Id == 0)
 	  {
 	    delete f;

src/loadGadget.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/loadGadget.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMO_LOAD_GADGET_HPP
 #define __COSMO_LOAD_GADGET_HPP
 
@@ -8,9 +43,8 @@
 
 namespace CosmoTool {
 
-  PurePositionData *loadGadgetPosition(const char *fname);
+  SimuData *loadGadgetMulti(const char *fname, int id, int flags, 
-
+			    int GadgetFormat = 1, SimuFilter filter = 0);
-  SimuData *loadGadgetMulti(const char *fname, int id, int flags, int GadgetFormat = 1);
  
   // Only single snapshot supported
   void writeGadget(const std::string& fname, SimuData *data, int GadgetFormat = 1);

src/loadRamses.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/loadRamses.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <sys/types.h>
 #include <regex.h>
 #include <cmath>
@@ -318,6 +353,9 @@ CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputI
   data->Omega_M = info.omega_m;
   data->Omega_Lambda = info.omega_lambda;
 
+  if (flags == 0)
+    return data;
+
   if (cpuid < 0)
     cpuid = 1;
   else cpuid++;
@@ -381,7 +419,7 @@ CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputI
 	}
       if (flags & NEED_GADGET_ID)
 	{
-	  data->Id = new int[nPar];
+	  data->Id = new long[nPar];
 	}
 
       for (int k = 0; k < 3; k++)

src/loadRamses.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/loadRamses.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef _LOAD_RAMSES_HPP
 #define _LOAD_RAMSES_HPP
 

src/loadSimu.hpp

@@ -1,6 +1,43 @@
+/*+
+This is CosmoTool (./src/loadSimu.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMOTOOLBOX_HPP
 #define __COSMOTOOLBOX_HPP
 
+#include <map>
+#include <string>
 
 namespace CosmoTool
 {
@@ -9,9 +46,24 @@ namespace CosmoTool
   static const int NEED_VELOCITY = 4;
   static const int NEED_TYPE = 8;
 
+  struct SimuParticle
+  {
+    float Pos[3];
+    float Vel[3];
+    int type;
+    int id;
+
+    bool flag_vel, flag_type, flag_id;
+  };
+
+  typedef bool (*SimuFilter)(const SimuParticle& p);
+
   class SimuData 
   {
   public:
+    typedef void (*FreeFunction)(void *);
+    typedef std::map<std::string, std::pair<void *, FreeFunction> > AttributeMap;
+
     float BoxSize;
     float time;
     float Hubble;
@@ -21,10 +73,13 @@ namespace CosmoTool
 
     long NumPart;
     long TotalNumPart;
-    int *Id;
+    long *Id;
     float *Pos[3];
     float *Vel[3];
     int *type;
+
+    AttributeMap attributes;
+    
   public:
     SimuData() : Id(0),NumPart(0),type(0)  { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; }
     ~SimuData() 
@@ -40,7 +95,37 @@ namespace CosmoTool
 	delete[] type;
       if (Id)
 	delete[] Id;
+
+      for (AttributeMap::iterator i = attributes.begin();
+	   i != attributes.end();
+	   ++i)
+	{
+	  if (i->second.second)
+	    i->second.second(i->second.first);
 	}
+    }    
+    
+    template<typename T>
+    T *as(const std::string& n) 
+    {
+      AttributeMap::iterator i = attributes.find(n);
+      if (i == attributes.end())
+	return 0;
+      
+      return reinterpret_cast<T *>(i->second.first);
+    }
+
+    void new_attribute(const std::string& n, void *p, FreeFunction free_func)
+    {
+      AttributeMap::iterator i = attributes.find(n);
+      if (i != attributes.end())
+	{
+	  if (i->second.second)
+	    i->second.second(i->second.first);
+	}
+      attributes[n] = std::make_pair(p, free_func);
+    }
+    
   };
 
 };

src/load_data.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/load_data.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <assert.h>
 #include <stdio.h>
 #include <stdlib.h>

src/load_data.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/load_data.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef _LOAD_GADGET_DATA_HPP
 #define _LOAD_GADGET_DATA_HPP
 

src/mach.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/mach.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMO_MACHINE_TEST_HPP
 #define __COSMO_MACHINE_TEST_HPP
 

src/miniargs.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/miniargs.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <cstdlib>
 #include <cstdio>
 #include <cstring>

src/miniargs.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/miniargs.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef _MAK_MINIARGS_HPP
 #define _MAK_MINIARGS_HPP
 

src/mykdtree.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/mykdtree.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __HV_KDTREE_HPP
 #define __HV_KDTREE_HPP
 
@@ -65,15 +100,16 @@ namespace CosmoTool {
   class RecursionMultipleInfo
   {
   public:
-    const typename KDDef<N,CType>::KDCoordinates& x;
+    typename KDDef<N,CType>::KDCoordinates x;
     BoundedQueue< KDCell<N,ValType,CType> *, typename KDDef<N,CType>::CoordType> queue;
     int traversed;
 
     RecursionMultipleInfo(const typename KDDef<N,CType>::KDCoordinates& rx,
 			  KDCell<N,ValType,CType> **cells,
 			  uint32_t numCells)
-      : x(rx), queue(cells, numCells, INFINITY),traversed(0)
+      : queue(cells, numCells, INFINITY),traversed(0)
     {
+      std::copy(rx, rx+N, x);
     }
   };
 
@@ -97,6 +133,12 @@ namespace CosmoTool {
     KDTree(Cell *cells, uint32_t Ncells);
     ~KDTree();
 
+    void setPeriodic(bool on, CoordType replicate)
+     {
+       periodic = on;
+       this->replicate = replicate;
+     }
+
     uint32_t getIntersection(const coords& x, CoordType r, 
 			     Cell **cells,
 			     uint32_t numCells)
@@ -150,6 +192,9 @@ namespace CosmoTool {
     Cell **sortingHelper;
     Cell *base_cell;
 
+    bool periodic;
+    CoordType replicate;
+
     Node *buildTree(Cell **cell0,
 		    uint32_t NumCells,
 		    uint32_t depth,

src/mykdtree.tcc

@@ -1,3 +1,4 @@
+#include "replicateGenerator.hpp"
 #include <cstring>
 #include <algorithm>
 #include <limits>
@@ -31,6 +32,7 @@ namespace CosmoTool {
   template<int N, typename ValType, typename CType, typename CellSplitter>
   KDTree<N,ValType,CType,CellSplitter>::KDTree(Cell *cells, uint32_t Ncells)
   {
+    periodic = false; 
 
     base_cell = cells;
     numNodes = Ncells;
@@ -93,6 +95,19 @@ namespace CosmoTool {
     info.distances = 0;
 
     recursiveIntersectionCells<false>(info, root, 0);
+    if (periodic)
+      {
+        ReplicateGenerator<float, N> r(x, replicate);
+
+        do
+          {
+            coords x_new;
+            r.getPosition(info.x);
+            recursiveIntersectionCells<false>(info, root, 0);
+          }
+        while (r.next());
+      }
+
     return info.currentRank;
   }
 
@@ -114,6 +129,18 @@ namespace CosmoTool {
     info.distances = distances;
 
     recursiveIntersectionCells<false>(info, root, 0);
+    if (periodic)
+      {
+        ReplicateGenerator<float, N> r(x, replicate);
+
+        do
+          {
+            coords x_new;
+            r.getPosition(info.x);
+            recursiveIntersectionCells<false>(info, root, 0);
+          }
+        while (r.next());
+      }
     return info.currentRank;
   }
 
@@ -131,6 +158,19 @@ namespace CosmoTool {
     info.distances = 0;
 
     recursiveIntersectionCells<true>(info, root, 0);
+    if (periodic)
+      {
+        ReplicateGenerator<float, N> r(x, replicate);
+
+        do
+          {
+            coords x_new;
+            r.getPosition(info.x);
+            recursiveIntersectionCells<true>(info, root, 0);
+          }
+        while (r.next());
+      }
+
     return info.currentRank;
   }
 
@@ -366,6 +406,18 @@ namespace CosmoTool {
    Cell *best = 0;
 
    recursiveNearest(root, 0, x, R2, best);
+   if (periodic)
+     {
+       ReplicateGenerator<float, N> r(x, replicate);
+
+       do
+         {
+           coords x_new;
+           r.getPosition(x_new);
+           recursiveNearest(root, 0, x_new, R2, best); 
+         }
+       while (r.next());
+     }
    
    return best;
  }
@@ -437,6 +489,18 @@ namespace CosmoTool {
      cells[i] = 0;
 
    recursiveMultipleNearest(info, root, 0);
+   if (periodic)
+    {
+       ReplicateGenerator<float, N> r(x, replicate);
+
+       do
+         {
+           coords x_new;
+           r.getPosition(info.x);
+           recursiveMultipleNearest(info, root, 0);
+         }
+       while (r.next());
+    }
 
    //   std::cout << "Traversed = " << info.traversed << std::endl;
  }
@@ -452,9 +516,19 @@ namespace CosmoTool {
      cells[i] = 0;
 
    recursiveMultipleNearest(info, root, 0);
-   memcpy(distances, info.queue.getPriorities(), sizeof(CoordType)*N2);
+   if (periodic)
+    {
+       ReplicateGenerator<float, N> r(x, replicate);
 
-   //   std::cout << "Traversed = " << info.traversed << std::endl;
+       do
+         {
+           coords x_new;
+           r.getPosition(info.x);
+           recursiveMultipleNearest(info, root, 0);
+         }
+       while (r.next());
+    }
+   memcpy(distances, info.queue.getPriorities(), sizeof(CoordType)*N2);
  }
 
 #ifdef __KD_TREE_SAVE_ON_DISK

src/newton.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/newton.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef _COSMOTOOL_NEWTON_HPP
 #define _COSMOTOOL_NEWTON_HPP
 

src/octTree.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/octTree.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include <iostream>
 #include <cmath>
 #include <cassert>

src/octTree.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/octTree.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMOTOOL_AMR_HPP
 #define __COSMOTOOL_AMR_HPP
 

src/pool.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/pool.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMO_POOL_HPP
 #define __COSMO_POOL_HPP
 

src/powerSpectrum.cpp

@@ -1,3 +1,39 @@
+/*+
+This is CosmoTool (./src/powerSpectrum.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
+#include <cassert>
 #include <vector>
 #include <algorithm>
 #include <iostream>
@@ -21,7 +57,7 @@ using namespace std;
 #define POWER_BDM 7
 #define POWER_TEST 8
 
-#define POWER_SPECTRUM HU_WIGGLES
+#define POWER_SPECTRUM POWER_EFSTATHIOU
 
 namespace Cosmology {
 
@@ -651,4 +687,17 @@ double computeCorrel2(double powNorm,  double topHatRad1, double topHatRad2)
 #endif
   }
 
+  double vvCorrection(double P_deltadelta, double k)
+  {
+    static const double alpha0 = -12480.5, alpha1 = 1.824, alpha2 = 2165.87, alpha3=1.796;
+    if (k > 0.3)
+      return 0;
+    double r =(alpha0*sqrt(P_deltadelta) + alpha1*P_deltadelta*P_deltadelta)/(alpha2 + alpha3*P_deltadelta);
+    assert(P_deltadelta > 0);
+
+    if (r < 0)
+      return 0;
+    return r;
+  }
+
 };

src/powerSpectrum.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/powerSpectrum.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef _POWERSPECTRUM_HPP
 #define _POWERSPECTRUM_HPP
 
@@ -30,6 +65,8 @@ namespace Cosmology {
   double computeVarianceZero(double powNorm);
   double computeCorrel(double powNorm,  double topHatRad1);
   double computeCorrel2(double powNorm,  double topHatRad1, double topHatRad2);
+
+  double vvCorrection(double P_deltadelta, double k);
 };
 
 #endif

src/replicateGenerator.hpp

@@ -0,0 +1,100 @@
+/*+
+This is CosmoTool (./src/replicateGenerator.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+#ifndef __REPLICATE_GENERATOR_HPP
+#define __REPLICATE_GENERATOR_HPP
+
+#include <algorithm>
+#include "algo.hpp"
+
+namespace CosmoTool
+{
+
+  template<typename Coord, int N>
+  class ReplicateGenerator
+  {
+  public:
+    typedef Coord Coords[N];
+    Coord replicate;
+
+    ReplicateGenerator(const Coords& x, Coord shift)
+    {
+      face = 0; 
+      replicate = shift;
+      numFaces = spower<N,long>(3);
+      std::copy(x, x+N, x_base);
+      if (!next())
+	abort();
+    }
+
+    bool next()
+    {
+      if (face == numFaces)
+        return false;
+
+      face++;
+
+      bool no_move = true;
+      int q_face = face;
+      for (int i = 0; i < N; i++)
+        {
+          int c_face;
+          c_face = q_face % 3;
+          q_face /= 3;
+          x_shifted[i] = x_base[i] + (c_face-1)*replicate;
+          no_move = no_move && (c_face == 1);
+        }
+      if (no_move)
+        return next();
+      return true;
+    }
+
+    const Coord *getPosition()
+    {
+      return x_shifted;
+    }
+
+    void getPosition(Coords& x_out)
+    {
+      std::copy(x_shifted, x_shifted+N, x_out);
+    }
+
+  private:
+    Coord x_shifted[N], x_base[N];
+    long face, numFaces;
+  };
+
+};
+
+#endif

src/sparseGrid.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/sparseGrid.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __VOID_SPARSE_GRID_HPP
 #define __VOID_SPARSE_GRID_HPP
 

src/sphSmooth.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/sphSmooth.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __COSMOTOOL_SPH_SMOOTH_HPP
 #define __COSMOTOOL_SPH_SMOOTH_HPP
 

src/yorick.hpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/yorick.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #ifndef __YORICK_HELPERS_HPP
 #define __YORICK_HELPERS_HPP
 

src/yorick_nc3.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/yorick_nc3.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include "ctool_netcdf_ver.hpp"
 #include "config.hpp"
 #ifdef NETCDFCPP4
@@ -151,7 +186,7 @@ namespace CosmoTool {
   ProgressiveOutput<T>::saveArrayProgressive(const std::string& fname, uint32_t *dimList,
 					     uint32_t rank)
   {
-    NcFile *f = new NcFile(fname.c_str(), NcFile::Replace);
+    NcFile *f = new NcFile(fname.c_str(), NcFile::Replace, 0, 0, NcFile::Netcdf4);
     
     assert(f->is_valid());
     
@@ -202,7 +237,7 @@ namespace CosmoTool {
   void saveArray(const std::string& fname,
 		 T *array, uint32_t *dimList, uint32_t rank)
   {
-    NcFile f(fname.c_str(), NcFile::Replace);
+    NcFile f(fname.c_str(), NcFile::Replace, 0, 0, NcFile::Netcdf4);
     
     assert(f.is_valid());
     

src/yorick_nc4.cpp

@@ -1,3 +1,38 @@
+/*+
+This is CosmoTool (./src/yorick_nc4.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013)
+
+guilhem.lavaux@gmail.com
+
+This software is a computer program whose purpose is to provide a toolbox for cosmological
+data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software.  You can  use, 
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info". 
+
+As a counterpart to the access to the source code and  rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty  and the software's author,  the holder of the
+economic rights,  and the successive licensors  have only  limited
+liability. 
+
+In this respect, the user's attention is drawn to the risks associated
+with loading,  using,  modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean  that it is complicated to manipulate,  and  that  also
+therefore means  that it is reserved for developers  and  experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or 
+data to be ensured and,  more generally, to use and operate it in the 
+same conditions as regards security. 
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
++*/
+
 #include "ctool_netcdf_ver.hpp"
 #include "config.hpp"
 #include <netcdf>