diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp index 9848597..3278876 100644 --- a/src/cosmopower.cpp +++ b/src/cosmopower.cpp @@ -7,16 +7,16 @@ This software is a computer program whose purpose is to provide a toolbox for co data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) This software is governed by the CeCILL license under French law and -abiding by the rules of distribution of free software. You can use, +abiding by the rules of distribution of free software. You can use, modify and/ or redistribute the software under the terms of the CeCILL license as circulated by CEA, CNRS and INRIA at the following URL -"http://www.cecill.info". +"http://www.cecill.info". As a counterpart to the access to the source code and rights to copy, modify and redistribute granted by the license, users are provided only with a limited warranty and the software's author, the holder of the economic rights, and the successive licensors have only limited -liability. +liability. In this respect, the user's attention is drawn to the risks associated with loading, using, modifying and/or developing or reproducing the @@ -25,9 +25,9 @@ that may mean that it is complicated to manipulate, and that also therefore means that it is reserved for developers and experienced professionals having in-depth computer knowledge. Users are therefore encouraged to load and test the software's suitability as regards their -requirements in conditions enabling the security of their systems and/or -data to be ensured and, more generally, to use and operate it in the -same conditions as regards security. +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. The fact that you are presently reading this means that you have had knowledge of the CeCILL license and that you accept its terms. @@ -54,7 +54,7 @@ using namespace CosmoTool; CosmoPower::CosmoPower() { eval = &CosmoPower::powerEfstathiou; - + n = 1.0; K0 = 1; V_LG_CMB = 627; @@ -71,7 +71,7 @@ CosmoPower::CosmoPower() Theta_27 = 2.728/2.7; ehu_params = 0; - + updateCosmology(); } @@ -86,7 +86,7 @@ CosmoPower::~CosmoPower() */ static double tophatFilter(double u) { - if (u != 0) + if (u != 0) return 3 / (u*u*u) * (sin(u) - u * cos(u)); else return 1; @@ -117,14 +117,14 @@ static double powG(double y) return y * (-6 * a + (2 + 3 * y) *log((a + 1)/(a - 1))); } - + double CosmoPower::powerEfstathiou(double k) { double a = 6.4/Gamma0; double b = 3/Gamma0; double c = 1.7/Gamma0; double nu = 1.13; - + double f = (a*k) + pow(b*k,1.5) + pow(c*k,2); // EFSTATHIOU ET AL. (1992) @@ -135,7 +135,7 @@ void CosmoPower::updateHuWigglesOriginal() { if (ehu_params == 0) ehu_params = new TF_Transfer(); - + ehu_params->TFset_parameters( (OMEGA_C+OMEGA_B)*h*h, OMEGA_B/(OMEGA_C+OMEGA_B), Theta_27*2.7); } @@ -161,7 +161,7 @@ void CosmoPower::updateHuWigglesConsts() 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, -f_b) * pow(a2, -pow(f_b, 3)); + double alpha_c = pow(a1, -f_b) * pow(a2, -pow(f_b, 3)); double b1_betac = 0.944 * 1/(1 + pow(458 * OmegaEff, -0.708)); double b2_betac = pow(0.395 * OmegaEff, -0.0266); @@ -170,13 +170,13 @@ void CosmoPower::updateHuWigglesConsts() 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 + f_b + (3 - 2 * f_b) * 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; @@ -194,24 +194,61 @@ double CosmoPower::powerHuWiggles(double k) double beta_c = ehu.beta_c; double alpha_b = ehu.alpha_b; double beta_b = ehu.beta_b; - double xx = k * s; + double xx = k * s; double s_tilde = s * pow(1 + pow(ehu.beta_node / (xx), 3), -1./3); - + double f = 1 / (1 + pow(xx / 5.4, 4)); 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(xx / 5.2, 2)) + + 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); + ) * j_0(k * s_tilde); - double T_k = OMEGA_B/OMEGA_0 * T_b + OMEGA_C/OMEGA_0 * T_c; + double T_k = OMEGA_B/OMEGA_0 * T_b + OMEGA_C/OMEGA_0 * T_c; return normPower * pow(k,n) * T_k * T_k; } +double CosmoPower::primordialPowerSpectrum(double k) +{ + //Primordial power spectrum, needed for PNG + + return normPower * pow(k,n); +} + +double CosmoPower::matterTransferFunctionHu(double k) +{ + // EISENSTEIN ET HU (1998) + // FULL POWER SPECTRUM WITH BARYONS AND WIGGLES + + double k_silk = ehu.k_silk; + double s = ehu.s; + double k_eq = ehu.k_eq; + double alpha_c = ehu.alpha_c; + double beta_c = ehu.beta_c; + double alpha_b = ehu.alpha_b; + double beta_b = ehu.beta_b; + double xx = k * s; + + double s_tilde = s * pow(1 + pow(ehu.beta_node / (xx), 3), -1./3); + + double f = 1 / (1 + pow(xx / 5.4, 4)); + 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(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; + + return T_k; +} + double CosmoPower::powerHuBaryons(double k) { double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75))); @@ -221,7 +258,7 @@ double CosmoPower::powerHuBaryons(double k) double L_0 = log(2 * M_E + 1.8 * q); double C_0 = 14.2 + 731 / (1 + 62.5 * q); double T0 = L_0 / (L_0 + C_0 * q * q); - + return normPower * pow(k,n) * T0 * T0; } @@ -350,7 +387,7 @@ double CosmoPower::eval_theta_theta(double k) return 0; double r =(alpha0*sqrt(P_deltadelta) + alpha1*P_deltadelta*P_deltadelta)/(alpha2 + alpha3*P_deltadelta); assert(P_deltadelta > 0); - + if (r < 0) return 0; return r; @@ -364,10 +401,10 @@ double CosmoPower::power(double k) double CosmoPower::powerHuWigglesOriginal(double k) { float tfb, tfc; - + ehu_params->TFfit_onek(k, &tfb, &tfc); - double T_k = OMEGA_B/OMEGA_0 * tfb + OMEGA_C/OMEGA_0 * tfc; + double T_k = OMEGA_B/OMEGA_0 * tfb + OMEGA_C/OMEGA_0 * tfc; return normPower * pow(k,n) * T_k * T_k; } @@ -383,6 +420,12 @@ void CosmoPower::setFunction(CosmoFunction f) updateHuWigglesConsts(); eval = &CosmoPower::powerHuWiggles; break; + case PRIMORDIAL_PS: + eval = &CosmoPower::primordialPowerSpectrum; + break; + case MATTER_TK: + eval = &CosmoPower::matterTransferFunctionHu; + break; case HU_WIGGLES_ORIGINAL: updateHuWigglesOriginal(); eval = &CosmoPower::powerHuWigglesOriginal; @@ -414,4 +457,3 @@ void CosmoPower::setNormalization(double A_K) { normPower = A_K;///power(0.002); } - diff --git a/src/cosmopower.hpp b/src/cosmopower.hpp index a04442b..ba0bffc 100644 --- a/src/cosmopower.hpp +++ b/src/cosmopower.hpp @@ -7,16 +7,16 @@ This software is a computer program whose purpose is to provide a toolbox for co data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) This software is governed by the CeCILL license under French law and -abiding by the rules of distribution of free software. You can use, +abiding by the rules of distribution of free software. You can use, modify and/ or redistribute the software under the terms of the CeCILL license as circulated by CEA, CNRS and INRIA at the following URL -"http://www.cecill.info". +"http://www.cecill.info". As a counterpart to the access to the source code and rights to copy, modify and redistribute granted by the license, users are provided only with a limited warranty and the software's author, the holder of the economic rights, and the successive licensors have only limited -liability. +liability. In this respect, the user's attention is drawn to the risks associated with loading, using, modifying and/or developing or reproducing the @@ -25,9 +25,9 @@ that may mean that it is complicated to manipulate, and that also therefore means that it is reserved for developers and experienced professionals having in-depth computer knowledge. Users are therefore encouraged to load and test the software's suitability as regards their -requirements in conditions enabling the security of their systems and/or -data to be ensured and, more generally, to use and operate it in the -same conditions as regards security. +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. The fact that you are presently reading this means that you have had knowledge of the CeCILL license and that you accept its terms. @@ -39,7 +39,7 @@ knowledge of the CeCILL license and that you accept its terms. namespace CosmoTool { struct TF_Transfer; - + class CosmoPower { public: @@ -57,15 +57,15 @@ namespace CosmoTool { double omega_B; double omega_C; double Theta_27; - - // DERIVED VARIABLES + + // DERIVED VARIABLES double OMEGA_0; double Omega; double beta; double OmegaEff; double Gamma0; double normPower; - + struct EHuParams { double k_silk; double s; @@ -74,13 +74,15 @@ namespace CosmoTool { double alpha_c, beta_c; double beta_node; }; - + EHuParams ehu; TF_Transfer *ehu_params; enum CosmoFunction { POWER_EFSTATHIOU, HU_WIGGLES, + PRIMORDIAL_PS, + MATTER_TK, HU_BARYON, OLD_POWERSPECTRUM, POWER_BARDEEN, @@ -101,7 +103,7 @@ namespace CosmoTool { void setNormalization(double A_K); void updateHuWigglesConsts(); void updateHuWigglesOriginal(); - + double eval_theta_theta(double k); double power(double k); @@ -111,6 +113,8 @@ namespace CosmoTool { double powerEfstathiou(double k); double powerHuWiggles(double k); + double primordialPowerSpectrum(double k); + double matterTransferFunctionHu(double k); double powerHuBaryons(double k); double powerOld(double k); double powerBardeen(double k);