From 72df7b6c6613b813edfb6db6f28f8fbdc50721e8 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem.lavaux@gmail.com>
Date: Thu, 13 Dec 2012 17:27:46 -0500
Subject: [PATCH] Fully implemented cosmopower

---
 src/cosmopower.cpp                     | 38 ++++++++++++++++++++++++--
 src/cosmopower.hpp                     |  4 +--
 src/fourier/details/euclidian_maps.hpp | 10 +++----
 3 files changed, 43 insertions(+), 9 deletions(-)

diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp
index 7f64080..e82c25e 100644
--- a/src/cosmopower.cpp
+++ b/src/cosmopower.cpp
@@ -197,7 +197,7 @@ double CosmoPower::integrandNormalize(double k)
   return power(k)*k*k*f*f;
 }
 
-void CosmoPower::normalize(double sigma8)
+void CosmoPower::normalize()
 {
   int Nsteps = 30000;
   double normVal = 0;
@@ -216,7 +216,7 @@ void CosmoPower::normalize(double sigma8)
 
   normVal /= (2*M_PI*M_PI);
 
-  normPower =  sigma8*sigma8/normVal;
+  normPower =  SIGMA8*SIGMA8/normVal;
 }
 
 void CosmoPower::updateCosmology()
@@ -248,3 +248,37 @@ 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();
+    }
+}
diff --git a/src/cosmopower.hpp b/src/cosmopower.hpp
index 83b3d8a..b700fc0 100644
--- a/src/cosmopower.hpp
+++ b/src/cosmopower.hpp
@@ -41,10 +41,10 @@ namespace CosmoTool {
 
     CosmoPower();
 
-    void setFunction();
+    void setFunction(CosmoFunction f);
 
     void updateCosmology();
-    void normalize(double sigma8);
+    void normalize();
 
  
     double eval_theta_theta(double k);
diff --git a/src/fourier/details/euclidian_maps.hpp b/src/fourier/details/euclidian_maps.hpp
index 6a650fe..a7c51b5 100644
--- a/src/fourier/details/euclidian_maps.hpp
+++ b/src/fourier/details/euclidian_maps.hpp
@@ -119,7 +119,7 @@ namespace CosmoTool
     }
 
     template<typename Array, typename Array2>
-    void get_Kvec(const Array& ik, Array2& kvec)
+    void get_Kvec(const Array& ik, Array2& kvec) const
     {
       const DimArray& dims = this->getDims();
       assert(ik.size() == dims.size());
@@ -137,7 +137,7 @@ namespace CosmoTool
     }
 
     template<typename Array2>
-    void get_Kvec(long p, Array2& kvec)
+    void get_Kvec(long p, Array2& kvec) const
     {
       const DimArray& dims = this->getDims();
       DimArray d(delta_k.size());
@@ -145,7 +145,7 @@ namespace CosmoTool
       get_Kvec(d, kvec);
     }
 
-    void get_IKvec(long p, DimArray& ikvec)
+    void get_IKvec(long p, DimArray& ikvec) const
     {
       const DimArray& dims = this->getDims();
       assert(dims.size()==ikvec.size());
@@ -158,7 +158,7 @@ namespace CosmoTool
 
 
     template<typename Array>
-    double get_K(const Array& ik)
+    double get_K(const Array& ik) const
     {
       const DimArray& dims = this->getDims();
       assert(ik.size() == dims.size());
@@ -177,7 +177,7 @@ namespace CosmoTool
       return std::sqrt(k2);
     }
 
-    double get_K(long p)
+    double get_K(long p) const
     {
       const DimArray& dims = this->getDims();
       DimArray d(delta_k.size());