From b7392204adc6ed7a378eeb10c060b29a8213cfda Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Tue, 24 Feb 2015 14:48:08 +0100
Subject: [PATCH] Added possibility to do partial normalization

---
 python/_cosmo_power.pyx | 6 +++---
 src/cosmopower.cpp      | 8 ++++++--
 src/cosmopower.hpp      | 2 +-
 3 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/python/_cosmo_power.pyx b/python/_cosmo_power.pyx
index 4e9e871..2e87c5c 100644
--- a/python/_cosmo_power.pyx
+++ b/python/_cosmo_power.pyx
@@ -44,7 +44,7 @@ cdef extern from "cosmopower.hpp" namespace "CosmoTool":
     void setFunction(CosmoFunction)
     void updateCosmology()
     void updatePhysicalCosmology()
-    void normalize()
+    void normalize(double)
     void setNormalization(double)
     double power(double)
  
@@ -75,7 +75,7 @@ cdef class CosmologyPower:
     
     self.power.updateCosmology()
     
-  def normalize(self,s8):
+  def normalize(self,s8,k_max=-1):
     """normalize(self, sigma8)
     
     Compute the normalization of the power spectrum using sigma8.
@@ -84,7 +84,7 @@ cdef class CosmologyPower:
       sigma8 (float): standard deviation of density field smoothed at 8 Mpc/h
     """
     self.power.SIGMA8 = s8
-    self.power.normalize()
+    self.power.normalize(k_max)
     
     
   def setFunction(self,funcname):
diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp
index 8e6b940..5104148 100644
--- a/src/cosmopower.cpp
+++ b/src/cosmopower.cpp
@@ -262,12 +262,16 @@ double CosmoPower::integrandNormalize(double x)
   return power(k)*k*k*f*f/(x*x);
 }
 
-void CosmoPower::normalize()
+void CosmoPower::normalize(double k_max)
 {
   double normVal = 0;
   double abserr;
   gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION);
   gsl_function f;
+  double x_min = 0;
+
+  if (k_max > 0)
+    x_min = 1/(1+k_max);
 
   f.function = gslPowSpecNorm;
   f.params = this;
@@ -282,7 +286,7 @@ void CosmoPower::normalize()
     }
 
  // 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_qag(&f, x_min, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr);
   gsl_integration_workspace_free(w);
 
   normVal /= (2*M_PI*M_PI);
diff --git a/src/cosmopower.hpp b/src/cosmopower.hpp
index bf2de61..f1ba45d 100644
--- a/src/cosmopower.hpp
+++ b/src/cosmopower.hpp
@@ -92,7 +92,7 @@ namespace CosmoTool {
 
     void updateCosmology();
     void updatePhysicalCosmology();
-    void normalize();
+    void normalize(double k_max = -1);
     void setNormalization(double A_K);
     void updateHuWigglesConsts();