Fixed bug in CosmoPower/ Eisenstein & Hu

python/_cosmo_cic.pyx

@@ -1,5 +1,6 @@
 from libcpp cimport bool
 from libcpp cimport string as cppstring
+from libcpp cimport vector as cppvector
 import numpy as np
 cimport numpy as np
 from cpython cimport PyObject, Py_INCREF
@@ -29,7 +30,7 @@ cdef extern from "cic.hpp" namespace "CosmoTool":
 @cython.cdivision(True)
 @cython.wraparound(False)
 def leanCic(float[:,:] particles, float L, int Resolution):
-  cdef CICParticles p
+  cdef cppvector.vector[CICParticles] *p
   cdef CICFilter *cic
   cdef np.uint64_t i
   cdef CICType *field
@@ -40,25 +41,34 @@ def leanCic(float[:,:] particles, float L, int Resolution):
   cdef np.uint64_t j
 
   cic = new CICFilter(Resolution, L)
+  print("Reset mesh")
   cic.resetMesh()
   
   if particles.shape[1] != 3:
     raise ValueError("Particles must be Nx3 array")
   
-  p.mass = 1
+  print("Inserting particles")
-  for i in xrange(particles.shape[0]):
+#  p = new cppvector.vector[CICParticles](particles.shape[0])
-    p.coords[0] = particles[i,0]
+#  for i in xrange(particles.shape[0]):
-    p.coords[1] = particles[i,1]
+#    *p[i].mass = 1
-    p.coords[2] = particles[i,2]
+#    *p[i].coords[0] = particles[i,0]
-    cic.putParticles(&p, 1)
+#    *p[i].coords[1] = particles[i,1]
+#    *p[i].coords[2] = particles[i,2]
+
+#  cic.putParticles(&p[0], particles.shape[0])
+  del p
     
+  print("Done")
   field = <CICType*>0
   dummyRes = 0
   cic.getDensityField(field, dummyRes)
+
+  print("Got to allocate a numpy %dx%dx%d" % (dummyRes, dummyRes,dummyRes))
   
   out_field = np.empty((dummyRes, dummyRes, dummyRes), dtype=np.float64)
   out_field0 = out_field.reshape(out_field.size)
   out_field_buf = out_field
+  print("Copy")
   for j in xrange(out_field_buf.size):
     out_field_buf[j] = field[j]
   

python/_cosmo_power.pyx

@@ -44,7 +44,7 @@ cdef extern from "cosmopower.hpp" namespace "CosmoTool":
     void setFunction(CosmoFunction)
     void updateCosmology()
     void updatePhysicalCosmology()
-    void normalize(double)
+    void normalize(double,double)
     void setNormalization(double)
     double power(double)
  
@@ -74,8 +74,11 @@ cdef class CosmologyPower:
     assert self.power.OMEGA_C > 0
     
     self.power.updateCosmology()
+
+  def setNormalization(self,A):
+    self.power.setNormalization(A)
     
-  def normalize(self,s8,k_max=-1):
+  def normalize(self,s8,k_min=-1,k_max=-1):
     """normalize(self, sigma8)
     
     Compute the normalization of the power spectrum using sigma8.
@@ -84,7 +87,7 @@ cdef class CosmologyPower:
       sigma8 (float): standard deviation of density field smoothed at 8 Mpc/h
     """
     self.power.SIGMA8 = s8
-    self.power.normalize(k_max)
+    self.power.normalize(k_min, k_max)
     
     
   def setFunction(self,funcname):

src/cic.cpp

@@ -186,6 +186,11 @@ void CICFilter::putParticles(CICParticles *particles, uint32_t N)
     for (long p = tUsedMin[t]; p < tUsedMax[t]; p++)
       densityGrid[p] += threadedDensity[t][p];
   }
+
+  for (int t = 0; t < threadUsed; t++)
+  {
+    delete[] threadedDensity;
+  }
 }
   
 void CICFilter::getDensityField(CICType*& field, uint32_t& res)

src/cosmopower.cpp

@@ -90,7 +90,8 @@ static double powC(double q, double alpha_c)
 
 static double T_tilde_0(double q, double alpha_c, double beta_c)
 {
-  double a = log(M_E + 1.8 * beta_c * q);
+  static const double c_E = 2.718282; //M_E;
+  double a = log(c_E + 1.8 * beta_c * q);
   return a / ( a  + powC(q, alpha_c) * q * q);
 }
 
@@ -123,6 +124,7 @@ double CosmoPower::powerEfstathiou(double k)
 
 void CosmoPower::updateHuWigglesConsts()
 {
+
   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)));
@@ -130,7 +132,7 @@ void CosmoPower::updateHuWigglesConsts()
 
   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 z_d = 1291 * pow(OmegaEff, 0.251) / (1 + 0.659 * pow(OmegaEff, 0.828)) * (1 + b1_zd * pow(OMEGA_B*h*h, b2_zd));
 
   double R_d = 31.5 * OMEGA_B * h * h * pow(Theta_27, -4) * 1e3 / z_d;
   double Req = 31.5 * OMEGA_B * h * h * pow(Theta_27, -4) * 1e3 / z_eq;
@@ -180,7 +182,10 @@ double CosmoPower::powerHuWiggles(double k)
   double q = k / (13.41 * k_eq);
   double T_c = f * T_tilde_0(q, 1, beta_c) + (1 - f) * T_tilde_0(q, alpha_c, beta_c);
 
-  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_b = ( 
+     T_tilde_0(q, 1, 1) / (1 + pow(xx / 5.2, 2)) + 
+     alpha_b / (1 + pow(beta_b / xx, 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;  
 
@@ -262,16 +267,18 @@ double CosmoPower::integrandNormalize(double x)
   return power(k)*k*k*f*f/(x*x);
 }
 
-void CosmoPower::normalize(double k_max)
+void CosmoPower::normalize(double k_min, 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;
+  double x_min = 0, x_max = 1;
 
   if (k_max > 0)
     x_min = 1/(1+k_max);
+  if (k_min > 0)
+    x_max = 1/(1+k_min);
 
   f.function = gslPowSpecNorm;
   f.params = this;
@@ -281,12 +288,12 @@ void CosmoPower::normalize(double k_max)
   ofstream ff("PP_k.txt");
   for (int i = 0; i < 100; i++)
     {
-       double k = pow(10.0, 4.0*i/100.-2);
+       double k = pow(10.0, 8.0*i/100.-4);
        ff << k << " " << power(k) << endl;
     }
 
  // gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr);
-  gsl_integration_qag(&f, x_min, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr);
+  gsl_integration_qag(&f, x_min, x_max, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr);
   gsl_integration_workspace_free(w);
 
   normVal /= (2*M_PI*M_PI);
@@ -371,5 +378,5 @@ void CosmoPower::setFunction(CosmoFunction f)
 
 void CosmoPower::setNormalization(double A_K)
 {
-  normPower = A_K/power(0.002);
+  normPower = A_K;///power(0.002);
 }

src/cosmopower.hpp

@@ -92,7 +92,7 @@ namespace CosmoTool {
 
     void updateCosmology();
     void updatePhysicalCosmology();
-    void normalize(double k_max = -1);
+    void normalize(double k_min = -1, double k_max = -1);
     void setNormalization(double A_K);
     void updateHuWigglesConsts();