Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BUG: Fixed wavenumber definition and amplitude scaling in Eisenstein & Hu model #445

Merged
merged 5 commits into from
Mar 5, 2021
Merged

BUG: Fixed wavenumber definition and amplitude scaling in Eisenstein & Hu model #445

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
35524f2
Fixed wavenumber definition in E&Hu model
jucordero Mar 3, 2021
5a01491
Added missing h^3 factor
jucordero Mar 3, 2021
760b5f2
Consistent use of wavenumber, precise c, corrected tests
jucordero Mar 4, 2021
bf3bef9
Flake8 errors
jucordero Mar 4, 2021
bef1f63
Apply suggestions from code review
rrjbca Mar 5, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 20 additions & 18 deletions skypy/power_spectrum/_eisenstein_hu.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@

from astropy.utils import isiterable
import numpy as np
from astropy import constants


__all__ = [
Expand Down Expand Up @@ -54,8 +55,8 @@ def transfer_with_wiggles(wavenumber, A_s, n_s, cosmology, kwmap=0.02):
>>> wavenumber = np.logspace(-3, 1, num=5, base=10.0)
>>> A_s, n_s = 2.1982e-09, 0.969453
>>> transfer_with_wiggles(wavenumber, A_s, n_s, Planck15, kwmap=0.02)
array([9.92144790e-01, 7.78548704e-01, 1.29998169e-01, 4.63863054e-03,
8.87918075e-05])
array([9.84445487e-01, 6.77428507e-01, 8.21122028e-02, 2.44250638e-03,
4.41412308e-05])

References
----------
Expand Down Expand Up @@ -84,7 +85,6 @@ def transfer_with_wiggles(wavenumber, A_s, n_s, cosmology, kwmap=0.02):
raise ValueError("Tcmb0 for input cosmology must be non-zero if"
"wiggles = True")

ak = wavenumber * h0
om0h2 = om0 * h0**2
ob0h2 = ob0 * h0**2
f_baryon = ob0 / om0
Expand Down Expand Up @@ -126,8 +126,8 @@ def transfer_with_wiggles(wavenumber, A_s, n_s, cosmology, kwmap=0.02):
beta_b = 0.5 + f_baryon + (3 - 2 * f_baryon) * np.sqrt((17.2 * om0h2)**2
+ 1.0)

q = ak / (13.41 * k_eq)
ks = ak * sound_horizon
q = wavenumber / (13.41 * k_eq)
ks = wavenumber * sound_horizon

T_c_ln_beta = np.log(np.e + 1.8 * beta_c * q)
T_c_ln_nobeta = np.log(np.e + 1.8 * q)
Expand All @@ -143,11 +143,12 @@ def f(a, b):
(1 - T_c_f) * f(T_c_ln_beta, T_c_C_alpha)

s_tilde = sound_horizon * (1 + (beta_node / ks)**3)**(-1 / 3)
ks_tilde = ak * s_tilde
ks_tilde = wavenumber * s_tilde

T_b_T0 = f(T_c_ln_nobeta, T_c_C_noalpha)
T_b_1 = T_b_T0 / (1 + (ks / 5.2)**2)
T_b_2 = alpha_b / (1 + (beta_b / ks)**3) * np.exp(-(ak / k_silk)**1.4)
T_b_2 = alpha_b / (1 + (beta_b / ks)**3) * np.exp(-(wavenumber
/ k_silk)**1.4)
rrjbca marked this conversation as resolved.
Show resolved Hide resolved
T_b = np.sinc(ks_tilde / np.pi) * (T_b_1 + T_b_2)

transfer = f_baryon * T_b + (1 - f_baryon) * T_c
Expand Down Expand Up @@ -190,8 +191,8 @@ def transfer_no_wiggles(wavenumber, A_s, n_s, cosmology):
>>> wavenumber = np.logspace(-3, 1, num=5, base=10.0)
>>> A_s, n_s = 2.1982e-09, 0.969453
>>> transfer_no_wiggles(wavenumber, A_s, n_s, Planck15)
array([9.91959695e-01, 7.84518347e-01, 1.32327555e-01, 4.60773671e-03,
8.78447096e-05])
array([9.84321214e-01, 6.83446969e-01, 8.15258290e-02, 2.42055117e-03,
4.36725897e-05])

References
----------
Expand All @@ -210,16 +211,16 @@ def transfer_no_wiggles(wavenumber, A_s, n_s, cosmology):
ob0 = cosmology.Ob0
h0 = cosmology.H0.value / 100
Tcmb0 = cosmology.Tcmb0.value
ak = wavenumber * h0
om0h2 = om0 * h0**2
f_baryon = ob0 / om0

alpha = 1 - 0.328 * np.log(431 * om0h2) * f_baryon + 0.38 * \
np.log(22.3 * om0h2) * f_baryon**2
sound = 44.5 * np.log(9.83 / om0h2) / \
np.sqrt(1 + 10 * (f_baryon * om0h2)**(0.75))
shape = om0h2 * (alpha + (1 - alpha) / (1 + (0.43 * ak * sound)**4))
aq = ak * (Tcmb0 / 2.7)**2 / shape
shape = om0h2 * (alpha + (1 - alpha) / (1 + (0.43 * wavenumber
* sound)**4))
rrjbca marked this conversation as resolved.
Show resolved Hide resolved
aq = wavenumber * (Tcmb0 / 2.7)**2 / shape
transfer = np.log(2 * np.e + 1.8 * aq) / \
(np.log(2 * np.e + 1.8 * aq) +
(14.2 + 731 / (1 + 62.5 * aq)) * aq * aq)
Expand Down Expand Up @@ -272,8 +273,8 @@ def eisenstein_hu(wavenumber, A_s, n_s, cosmology, kwmap=0.02, wiggle=True):
>>> wavenumber = np.logspace(-3, 1, num=5, base=10.0)
>>> A_s, n_s = 2.1982e-09, 0.969453
>>> eisenstein_hu(wavenumber, A_s, n_s, Planck15, kwmap=0.02, wiggle=True)
array([6.47460158e+03, 3.71610099e+04, 9.65702614e+03, 1.14604456e+02,
3.91399918e-01])
array([2.99126326e+04, 1.32023496e+05, 1.80797616e+04, 1.49108261e+02,
4.53912529e-01])

References
----------
Expand All @@ -282,16 +283,17 @@ def eisenstein_hu(wavenumber, A_s, n_s, cosmology, kwmap=0.02, wiggle=True):
.. [3] Komatsu et al., ApJS, 180, 330 (2009)
"""
om0 = cosmology.Om0
h0 = cosmology.H0.value / 100
H0_c = (cosmology.H0 / constants.c).to_value('Mpc-1')

if wiggle:
transfer = transfer_with_wiggles(wavenumber, A_s, n_s, cosmology,
kwmap)
else:
transfer = transfer_no_wiggles(wavenumber, A_s, n_s, cosmology)

# Eq [74] in [3]
power_spectrum = A_s * (2 * wavenumber**2 * 2998**2 / 5 / om0)**2 * \
transfer**2 * (wavenumber * h0 / kwmap)**(n_s - 1) * 2 * \
np.pi**2 / wavenumber**3
power_spectrum = A_s * (2 * (wavenumber / H0_c)**2 / 5 / om0)**2 * \
transfer**2 * (wavenumber / kwmap)**(n_s - 1) * 2 * \
np.pi**2 / (wavenumber)**3

return power_spectrum
12 changes: 6 additions & 6 deletions skypy/power_spectrum/tests/test_eisenstein_hu.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,13 +37,13 @@ def test_eisenstein_hu():
wiggle=True)
pk_eisensteinhu_nw = eisenstein_hu(wavenumber, A_s, n_s, cosmology, kwmap,
wiggle=False)
pk_cosmosis_w = np.array([6.47460158e+03, 3.71610099e+04, 9.65702614e+03,
1.14604456e+02, 3.91399918e-01])
pk_cosmosis_nw = np.array([6.47218600e+03, 3.77330704e+04, 1.00062077e+04,
1.13082980e+02, 3.83094714e-01])
pk_pre_w = np.array([2.99126326e+04, 1.32023496e+05, 1.80797616e+04,
1.49108261e+02, 4.53912529e-01])
pk_pre_nw = np.array([2.99050810e+04, 1.34379783e+05, 1.78224637e+04,
1.46439700e+02, 4.44325443e-01])

assert np.allclose(pk_eisensteinhu_w, pk_cosmosis_w)
assert np.allclose(pk_eisensteinhu_nw, pk_cosmosis_nw)
assert np.allclose(pk_eisensteinhu_w, pk_pre_w)
assert np.allclose(pk_eisensteinhu_nw, pk_pre_nw)

# Test for failure when wavenumber <= 0
negative_wavenumber_scalar = 0
Expand Down