From 031284bdbc3802293e29b19971bc7ce0be84922d Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Thu, 15 Jan 2009 16:42:55 -0600 Subject: [PATCH] Included standard power spectra --- src/powerSpectrum.cpp | 654 ++++++++++++++++++++++++++++++++++++++++++ src/powerSpectrum.hpp | 34 +++ 2 files changed, 688 insertions(+) create mode 100644 src/powerSpectrum.cpp create mode 100644 src/powerSpectrum.hpp diff --git a/src/powerSpectrum.cpp b/src/powerSpectrum.cpp new file mode 100644 index 0000000..b73d019 --- /dev/null +++ b/src/powerSpectrum.cpp @@ -0,0 +1,654 @@ +#include +#include +#include +#include +#include +#include +#include "powerSpectrum.hpp" + +using namespace std; + +#define USE_GSL +#define TOLERANCE 1e-6 +#define NUM_ITERATION 8000 + +#define POWER_EFSTATHIOU 1 +#define HU_WIGGLES 2 +#define HU_BARYON 3 +#define OLD_POWERSPECTRUM 4 +#define POWER_BARDEEN 5 +#define POWER_SUGIYAMA 6 +#define POWER_BDM 7 +#define POWER_TEST 8 + +#define POWER_SPECTRUM HU_BARYON + +namespace Cosmology { + +double n = 1.0; +double K0 = 1; +double V0 = 627; + + double CMB_VECTOR[3] = { + 56.759, + -540.02, + 313.50 + }; + + +// WMAP5 +double h = 0.719; +double SIGMA8 = 0.77; +double OMEGA_B = 0.043969; +double OMEGA_C = 0.21259; + + +// WMAP5-modification +//double h = 0.719; +//double SIGMA8 = 0.77; +//double OMEGA_B = 0; +//double OMEGA_C = 0.21259+0.043969; + +// LCDM STRAUSS ? +//double h = 0.67; +//double SIGMA8 = 0.67; +//double OMEGA_B = 0; +//double OMEGA_C = 0.30; + +// SCDM STRAUSS +//double h = 0.5; +//double SIGMA8= 1.05; +//double OMEGA_B = 0; +//double OMEGA_C = 1; + +// Sugiyama test +//double h = 0.5; +//double SIGMA8= 0.5;//1.05; +//double OMEGA_B = 0.0125*4; +//double OMEGA_C = 0.1-OMEGA_B; + +// HU TEST +//double h = 0.5; +//double SIGMA8 = 0.5; +//double OMEGA_B = 0.09; +//double OMEGA_C = 0.21; + +// HDM STRAUSS +//double h = 0.5; +//double SIGMA8 = 0.86; +//double OMEGA_B = 0; +//double OMEGA_C = 1; + +// FOR "BEST FIT" +//double h = 0.82; +//double SIGMA8 = 0.76; +//double OMEGA_B = 0.043969; +//double OMEGA_C = 0.15259; + + +// FOR JUSZKIEWICZ CHECKING (CDM) ! WARNING ! He smoothes +// with a gaussian filter the density field, i.e. one has +// to multiply P(k) by exp(-k^2 R^2) with R the radius +// of the filter. Dammit ! +//double h = 0.5; +//double SIGMA8=1/2.5; +//double OMEGA_B=0.; +//double OMEGA_C=1; +//#define JUSZKIEWICZ_PATCH +//#define RJUSZ 6.0 + +// (BDM) +//double h = 0.5; +//double SIGMA8=1; +//double OMEGA_B=0.0; +//double OMEGA_C=0.4; + + +// FOR HU CHECKING +//double h = 0.5; +//double SIGMA8= 1; +//double OMEGA_B=0.09; +//double OMEGA_C=0.21; + + +double OMEGA_0 = OMEGA_B+OMEGA_C; +double Omega = OMEGA_0; +double Theta_27 = 2.728 / 2.7; +double beta = pow(OMEGA_0, 5./9); +double OmegaEff = OMEGA_0 * h * h; +double Gamma0 = OMEGA_0 * h * h; + + +/* + * This is \hat{tophat} + */ +double tophatFilter(double u) +{ + if (u != 0) + return 3 / (u*u*u) * (sin(u) - u * cos(u)); + else + return 1; +} + +/* + * This is \tilde{w} + */ +double externalFilter(double u) +{ + if (u != 0) + return 1 - sin(u)/u; + return 0.; +} + +double powC(double q, double alpha_c) +{ + return 14.2 / alpha_c + 386 / (1 + 69.9 * pow(q, 1.08)); +} + +double T_tilde_0(double q, double alpha_c, double beta_c) +{ + double a = log(M_E + 1.8 * beta_c * q); + return a / ( a + powC(q, alpha_c) * q * q); +} + +double j_0(double x) +{ + if (x == 0) + return 1.0; + return sin(x)/x; +} + +double powG(double y) +{ + double a = sqrt(1 + y); + + return y * (-6 * a + (2 + 3 * y) *log((a + 1)/(a - 1))); +} + +/* + * This function returns the power spectrum evaluated at k (in Mpc, not in Mpc/h). + */ +double powerSpectrum(double k, double normPower) +{ +#if POWER_SPECTRUM == POWER_EFSTATHIOU + 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) + return normPower * pow(k,n) * pow(1+pow(f,nu),(-2/nu)); +#endif + + // EISENSTEIN ET HU (1998) + // FULL POWER SPECTRUM WITH BARYONS AND WIGGLES + +#if POWER_SPECTRUM == HU_WIGGLES + // 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 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 f = 1 / (1 + pow(k * s / 5.4, 4)); + double k_eq = 7.46e-2 * OmegaEff * pow(Theta_27, -2); + 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 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; + + return normPower * pow(k,n) * T_k * T_k; +#endif + + // EISENSTEIN ET AL. (2008), SHAPED POWER SPECTRUM WITH BARYON, WITHOUT WIGGLES +#if POWER_SPECTRUM == HU_BARYON + double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75))); + double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2); + double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4))); + double q = k/(h*GammaEff) * pow(Theta_27, 2); + 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; +#endif + +#if POWER_SPECTRUM == OLD_POWERSPECTRUM + // OLD FUNCTION: + static const double l = 1/(Omega * h*h); + static const double alpha = 1.7 * l, beta = 9.0 * pow(l, 1.5), gamma = l*l; + return normPower * pow(k,n) * pow(1 + alpha * k + beta * pow(k,1.5) + gamma *k*k,-2); +#endif + +#if POWER_SPECTRUM == POWER_SUGIYAMA + double q = k * Theta_27*Theta_27 / (OmegaEff * exp(-OMEGA_B - sqrt(h/0.5)*OMEGA_B/OMEGA_0)); + double L0 = log(2*M_E + 1.8 * q); + double C0 = 14.2 + 731 / (1 + 62.5 * q); + double T_k = L0 / (L0 + C0 * q*q); + + // double poly = 1 + 3.89 * q + pow(16.1*q,2) + pow(5.46*q,3) + pow(6.71*q,4); + // double T_k = log(1+2.34*q)/(2.34*q) * pow(poly,-0.25); + + return normPower * pow(k,n) * T_k * T_k; + +#endif + +#if POWER_SPECTRUM == POWER_BARDEEN + double q = k / (OmegaEff); + + double poly = 1 + 3.89 * q + pow(16.1*q,2) + pow(5.46*q,3) + pow(6.71*q,4); + double T_k = log(1+2.34*q)/(2.34*q) * pow(poly,-0.25); + + return normPower * pow(k,n) * T_k * T_k; + +#endif + +#if POWER_SPECTRUM == POWER_BDM + + k /= h*h; + double alpha1 = 190; + double Gmu = 4.6; + double alpha2 = 11.5; + double alpha3 = 11; + double alpha4 = 12.55; + double alpha5 = 0.0004; + return normPower*k*alpha1*alpha1*Gmu*Gmu/(1+(alpha2*k)+pow(alpha3*k,2)+pow(alpha4*k,3))*pow(1+pow(alpha5/k,2), -2); +#endif + +#if POWER_SPECTRUM == POWER_TEST + return 1/(1+k*k); +#endif +} + +/* + * This function computes the normalization of the power spectrum. It requests + * a sigma8 (density fluctuations within 8 Mpc/h) + */ +double gslPowSpecNorm(double k, void *params) +{ + double f = tophatFilter(k*8.0/h); + return powerSpectrum(k, 1.0)*k*k*f*f; +} + +double computePowSpecNorm(double sigma8) +{ + int Nsteps = 30000; + double normVal = 0; + +#ifndef USE_GSL + for (int i = 1; i <= Nsteps; i++) + { + double t = i * 1.0/(Nsteps+1); + // Change of variable ! + double k = (1-t)/t * K0; + + // The filter + double filter_val = tophatFilter(k*8.0/h); + // The powerspectrum + double powVal = powerSpectrum(k, 1.0); + // Multiply by the tophat filter + powVal *= filter_val*filter_val; + + powVal *= k*k; + + // Account for change of variable + powVal /= (t*t); + + // Integrate ! + normVal += powVal; + } + + normVal /= 2*M_PI*M_PI; + + // The dt element + normVal *= 1.0/(Nsteps+1) * K0; +#else + double abserr; + gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION); + gsl_function f; + + f.function = gslPowSpecNorm; + f.params = 0; + + gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr); + + gsl_integration_workspace_free(w); + + normVal /= (2*M_PI*M_PI); +#endif + + return sigma8*sigma8/normVal; +} + +/* + * This function computes the variance of the Local Group velocity components + * for a survey which depth is topHatRad1 (in Mpc/h). This variance should + * be multiplied by (H \beta)^2 to be equal to a velocity^2. + */ +double gslVariance(double k, void *params) +{ + double R1 = *(double *)params; + double f = externalFilter(k * R1 / h); + + double a = f*f; + +#ifdef JUSZKIEWICZ_PATCH + a *= exp(-k*k*(RJUSZ*RJUSZ/(h*h))); +#endif + a *= powerSpectrum(k, 1.0); + + return a; +} + +double computeVariance(double powNorm, double topHatRad1) +{ + int Nsteps = 100000; + double varVal = 0; + +#ifndef USE_GSL + + for (int i = 1; i <= Nsteps; i++) + { + double t = i * 1.0/(Nsteps+1); + // Change of variable ! + double k = (1-t)/t * K0; + double powVal = powerSpectrum(k, powNorm); + + double filter1Val = externalFilter(k*topHatRad1/h); + +#ifdef JUSZKIEWICZ_PATCH + powVal *= exp(-k*k*(RJUSZ*RJUSZ/(h*h))); +#endif + powVal *= filter1Val*filter1Val; + + powVal /= (t*t); + + varVal += powVal; + } + + varVal *= 1.0/(Nsteps) * K0; + varVal *= 1.0/(6*M_PI*M_PI); + +#else + + double abserr; + gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION); + gsl_function f; + + f.function = gslVariance; + f.params = &topHatRad1; + + gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &varVal, &abserr); + + gsl_integration_workspace_free(w); + + varVal *= powNorm/(6*M_PI*M_PI); + +#endif + + return varVal; +} + +/* + * This function computes the same quantity as computeVariance but + * for a survey infinitely deep. + */ +double gslVarianceZero(double k, void *params) +{ + double a = 1.0; + +#ifdef JUSZKIEWICZ_PATCH + a *= exp(-k*k*(RJUSZ*RJUSZ/(h*h))); +#endif + a *= powerSpectrum(k, 1.0); + + return a; +} + +double computeVarianceZero(double powNorm) +{ + int Nsteps = 100000; + double varVal = 0; + +#ifndef USE_GSL + + for (int i = 1; i <= Nsteps; i++) + { + double t = i * 1.0/(Nsteps+1); + // Change of variable ! + double k = (1-t)/t * K0; + double powVal = powerSpectrum(k, powNorm); + +#ifdef JUSZKIEWICZ_PATCH + powVal *= exp(-k*k*(RJUSZ*RJUSZ/h*h)); +#endif + + powVal /= (t*t); + + varVal += powVal; + } + + varVal *= 1.0/(Nsteps+1) * K0; + varVal *= 1.0/(6*M_PI*M_PI); + +#else + + double abserr; + gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION); + gsl_function f; + + f.function = gslVarianceZero; + f.params = 0; + + gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &varVal, &abserr); + + gsl_integration_workspace_free(w); + + varVal *= powNorm/(6*M_PI*M_PI); +#endif + + return varVal; +} + +/* + * This function computes the correlation between the infinitely deep + * velocity of the Local Group and the one estimated from a survey + * which depth is topHatRad1. + * This corresponds to \gamma. + * This quantity must be multiplied by H \beta to be equal to a velocity^2. + */ +double gslCorrel(double k, void *params) +{ + double R1 = ((double *)params)[0]; + double a = externalFilter(k * R1 / h);// * externalFilter(k * R2 / h); + +#ifdef JUSZKIEWICZ_PATCH + a *= exp(-k*k*(RJUSZ*RJUSZ/(h*h))); +#endif + a *= powerSpectrum(k, 1.0); + + return a; +} + +double gslCorrelBis(double t, void *params) +{ + double k = (1-t)/t; + double v = gslCorrel(k, params); + + return v/(t*t); +} + +double gslCorrel2(double k, void *params) +{ + double R1 = ((double *)params)[0]; + double R2 = ((double *)params)[1]; + double a = externalFilter(k * R1 / h) * externalFilter(k * R2 / h); + +#ifdef JUSZKIEWICZ_PATCH + a *= exp(-k*k*(RJUSZ*RJUSZ/(h*h))); +#endif + a *= powerSpectrum(k, 1.0) ; + + return a; +} + +double gslCorrel2bis(double t, void *params) +{ + double k = (1-t)/t; + double v = gslCorrel2(k, params); + + return v/(t*t); +} + + +double computeCorrel(double powNorm, double topHatRad1) +{ + int Nsteps = 100000; + double varVal = 0; + +#ifndef USE_GSL + + for (int i = 1; i <= Nsteps; i++) + { + double t = i * 1.0/(Nsteps+1); + // Change of variable ! + double k = (1-t)/t * K0; + double powVal = powerSpectrum(k, powNorm); + + double filter1Val = externalFilter(k*topHatRad1/h); + +#ifdef JUSZKIEWICZ_PATCH + powVal*=exp(-k*k*(RJUSZ*RJUSZ/h*h)); +#endif + + powVal *= filter1Val; + + powVal /= (t*t); + + varVal += powVal; + } + + varVal *= 1.0/(Nsteps) * K0; + varVal *= 1.0/(6*M_PI*M_PI); +#else + + double abserr; + gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION); + gsl_function f; + + f.params = &topHatRad1; + +#if 1 + f.function = gslCorrelBis; + gsl_integration_qag (&f, 0, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &varVal, &abserr); +#else + f.function = gslCorrel; + gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &varVal, &abserr); +#endif + + gsl_integration_workspace_free(w); + + varVal *= powNorm/(6*M_PI*M_PI); + +#endif + return varVal; +} + +/* + * This function computes the correlation between the infinitely deep + * velocity of the Local Group and the one estimated from a survey + * which depth is topHatRad1. + * This corresponds to \gamma. + * This quantity must be multiplied by H \beta to be equal to a velocity^2. + */ +double computeCorrel2(double powNorm, double topHatRad1, double topHatRad2) +{ + int Nsteps = 100000; + double varVal = 0; +#ifndef USE_GSL + + for (int i = 1; i <= Nsteps; i++) + { + double t = i * 1.0/(Nsteps+1); + // Change of variable ! + double k = (1-t)/t * K0; + double powVal = powerSpectrum(k, powNorm); + + double filter1Val = externalFilter(k*topHatRad1/h); + double filter2Val = externalFilter(k*topHatRad2/h); + + powVal *= filter1Val * filter2Val; + + powVal /= (t*t); + + varVal += powVal; + } + + varVal *= 1.0/(Nsteps) * K0; + varVal *= 1.0/(6*M_PI*M_PI); + +#else + + double abserr; + gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION); + gsl_function f; + double rads[] = {topHatRad1, topHatRad2 }; + + f.params = rads; + +#if 1 + f.function = gslCorrel2bis; + gsl_integration_qag (&f, 0, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &varVal, &abserr); +#else + f.function = gslCorrel2; + gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &varVal, &abserr); +#endif + + gsl_integration_workspace_free(w); + + varVal *= powNorm/(6*M_PI*M_PI); + +#endif + + return varVal; +} + + + void updateCosmology() + { + OMEGA_0 = OMEGA_B+OMEGA_C; + Omega = OMEGA_0; + Theta_27 = 2.728 / 2.7; + beta = pow(OMEGA_0, 5./9); + OmegaEff = OMEGA_0 * h * h; + Gamma0 = OMEGA_0 * h * h; + +#if 0 + cout << "Cosmology is :" << endl + << " O0=" << OMEGA_0 << " Theta=" << Theta_27 << " beta=" << beta << " h=" << h << " G0=" << Gamma0 << endl + << " OmegaB=" << OMEGA_B << " Omega_C=" << OMEGA_C << endl; +#endif + } + +}; diff --git a/src/powerSpectrum.hpp b/src/powerSpectrum.hpp new file mode 100644 index 0000000..58b3c30 --- /dev/null +++ b/src/powerSpectrum.hpp @@ -0,0 +1,34 @@ +#ifndef _POWERSPECTRUM_HPP +#define _POWERSPECTRUM_HPP + +namespace Cosmology { + + extern double n; + extern double K0; + extern double V0; + + extern double CMB_VECTOR[3]; + + // WMAP5 + extern double h; + extern double SIGMA8; + extern double OMEGA_B; + extern double OMEGA_C; + + extern double OMEGA_0; + extern double Omega; + extern double Theta_27; + extern double beta; + extern double OmegaEff; + extern double Gamma0; + + void updateCosmology(); + double powerSpectrum(double k, double normPower); + double computePowSpecNorm(double sigma8); + double computeVariance(double powNorm, double topHatRad1); + double computeVarianceZero(double powNorm); + double computeCorrel(double powNorm, double topHatRad1); + double computeCorrel2(double powNorm, double topHatRad1, double topHatRad2); +}; + +#endif