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();