Preinit Hu-Wiggles constants. Fixed sound horizon.

src/cosmopower.cpp

@@ -121,36 +121,64 @@ double CosmoPower::powerEfstathiou(double k)
   return normPower * pow(k,n) * pow(1+pow(f,nu),(-2/nu));
 }
 
+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)));
+  double k_eq = 7.46e-2 * OmegaEff * pow(Theta_27, -2);
+
+  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 Req = 31.5 * OMEGA_B * h * h * pow(Theta_27, -4) * 1e3 / z_eq;
+
+  double s = 2./(3.*k_eq) * sqrt(6/Req) * log((sqrt(1 + R_d) + sqrt(R_d + Req))/(1 + sqrt(Req)));
+
+  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 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 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);
+  
+  ehu.k_silk = k_silk;
+  ehu.s = s;
+  ehu.k_eq = k_eq;
+  ehu.alpha_c = alpha_c;
+  ehu.beta_c = beta_c;
+  
+  ehu.alpha_b = alpha_b;
+  ehu.beta_b = beta_b;
+  ehu.beta_node = beta_node;
+}
+
 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 k_silk = ehu.k_silk;
-  double z_eq = 2.50e4 * OmegaEff * pow(Theta_27, -4);
+  double s = ehu.s;
-  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double k_eq = ehu.k_eq;
-  double f = 1 / (1 + pow(k * s / 5.4, 4));
+  double alpha_c = ehu.alpha_c;
-  double k_eq = 7.46e-2 * OmegaEff * pow(Theta_27, -2);
+  double beta_c = ehu.beta_c;
-  double a1 = pow(46.9 * OmegaEff, 0.670) * (1 + pow(32.1 * OmegaEff, -0.532));
+  double alpha_b = ehu.alpha_b;
-  double a2 = pow(12.0 * OmegaEff, 0.424) * (1 + pow(45.0 * OmegaEff, -0.582));
+  double beta_b = ehu.beta_b;
-  double alpha_c = pow(a1, -OMEGA_B/ OMEGA_0) * pow(a2, -pow(OMEGA_B / OMEGA_0, 3));  
+  double s_tilde = s * pow(1 + pow(ehu.beta_node / (k * s), 3), -1./3);
-
+  
+  double xx = k * s;  
+  double f = 1 / (1 + pow(xx / 5.4, 4));
   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;  
@@ -310,6 +338,7 @@ void CosmoPower::setFunction(CosmoFunction f)
       eval = &CosmoPower::powerEfstathiou;
       break;
     case HU_WIGGLES:
+      updateHuWigglesConsts();
       eval = &CosmoPower::powerHuWiggles;
       break;
     case HU_BARYON:

src/cosmopower.hpp

@@ -63,18 +63,28 @@ namespace CosmoTool {
     double OmegaEff;
     double Gamma0;
     double normPower;
+    
+    struct EHuParams {
+      double k_silk;
+      double s;
+      double k_eq;
+      double alpha_b, beta_b;
+      double alpha_c, beta_c;
+      double beta_node;
+    };
+    
+    EHuParams ehu;
 
-    enum CosmoFunction
+    enum CosmoFunction {
-      {
+      POWER_EFSTATHIOU,
-	POWER_EFSTATHIOU,
+      HU_WIGGLES,
-	HU_WIGGLES,
+      HU_BARYON,
-	HU_BARYON,
+      OLD_POWERSPECTRUM,
-	OLD_POWERSPECTRUM,
+      POWER_BARDEEN,
-	POWER_BARDEEN,
+      POWER_SUGIYAMA,
-	POWER_SUGIYAMA,
+      POWER_BDM,
-	POWER_BDM,
+      POWER_TEST
-	POWER_TEST
+    };
-      };
 
     CosmoPower();
 
@@ -84,7 +94,7 @@ namespace CosmoTool {
     void updatePhysicalCosmology();
     void normalize();
     void setNormalization(double A_K);
-
+    void updateHuWigglesConsts();
  
     double eval_theta_theta(double k);
     double power(double k);