Update constants, nnu=3.044, Mpc, Gyr, BBN model default to PRIMAT2021

.travis.yml

@@ -34,6 +34,15 @@ jobs:
         - PYDIST="ANACONDA"
         - FORTRAN="test"
       python: "3.10"
+    - name: "Ubuntu 22.04 + Python 3.11"
+      dist: jammy
+      addons:
+        apt:
+          packages:
+            - gfortran
+      env:
+        - CHANNEL="defaults"
+      python: "3.11"
     - name: "GCC trunk cosmobox"
       language: shell
       services:

camb/__init__.py

@@ -7,7 +7,7 @@
 __author__ = "Antony Lewis"
 __contact__ = "antony at cosmologist dot info"
 __url__ = "https://camb.readthedocs.io"
-__version__ = "1.3.8"
+__version__ = "1.4.0"
 
 from . import baseconfig
 

camb/bbn.py

@@ -2,41 +2,23 @@
 # Fitting formula for Parthenelope via Julien Lesgourgues Dec 2014
 # General interpolation tables, with defaults from Parthenope May 2017 (thanks Ofelia Pesanti)
 # Use PRIMAT_Yp_DH_Error.dat table for latest from the PRIMAT code (arXiv: 1801.08023, thanks Cyril Pitrou)
+# from 1.4.0, delta_nnu values are relative to 3.044, not 3.046
 
 import numpy as np
 import os
 from scipy.interpolate import RectBivariateSpline
+from .constants import m_H, m_He4, default_nnu
 
-# Various useful constants
-hbar = 1.05457e-34
-c = 299792458.
-kB = 1.380650e-23
-MeV = 1.6021764e-13
-eV = MeV / 1e6
-G = 6.6738e-11
-TCMB = 2.7255
-mP = np.sqrt(hbar * c / G)
-Mpc = 3.0856776e22
-
-m_proton = 1.672621637e-27
-m_H = 1.673575e-27
-not4 = 3.9715
-m_He = m_H * not4
-
-zeta3 = 1.202056903
-
-n_photon = (kB * TCMB / hbar / c) ** 3 * zeta3 * 2 / np.pi ** 2
-omegafac = (1e5 / Mpc) ** 2 / (8 * np.pi * G) * 3
-
-default_interpolation_table = 'PArthENoPE_880.2_standard.dat'
+# previously (< 1.4.0) default_interpolation_table = 'PArthENoPE_880.2_standard.dat'
+default_interpolation_table = 'PRIMAT_Yp_DH_ErrorMC_2021.dat'
 
 
 def yhe_to_ypBBN(Yp):
-    return -4 * m_H * Yp / (Yp * m_He - 4 * Yp * m_H - m_He)
+    return -4 * m_H * Yp / (Yp * m_He4 - 4 * Yp * m_H - m_He4)
 
 
 def ypBBN_to_yhe(YBBN):
-    return -YBBN * m_He / (-YBBN * m_He + 4 * YBBN * m_H - 4 * m_H)
+    return -YBBN * m_He4 / (-YBBN * m_He4 + 4 * YBBN * m_H - 4 * m_H)
 
 
 class BBNIterpolator(RectBivariateSpline):
@@ -53,7 +35,7 @@ def Y_p(self, ombh2, delta_neff=0.):
         Get BBN helium nucleon fraction. Must be implemented by extensions.
 
         :param ombh2: :math:`\Omega_b h^2`
-        :param delta_neff:  additional N_eff relative to standard value (of 3.046)
+        :param delta_neff:  additional N_eff relative to standard value (of 3.044)
         :return:  Y_p helium nucleon fraction predicted by BBN
         """
         raise Exception('Not implemented')
@@ -63,7 +45,7 @@ def Y_He(self, ombh2, delta_neff=0.):
         Get BBN helium mass fraction for CMB code.
 
         :param ombh2: :math:`\Omega_b h^2`
-        :param delta_neff:  additional N_eff relative to standard value (of 3.046)
+        :param delta_neff:  additional N_eff relative to standard value (of 3.044)
         :return: Y_He helium mass fraction predicted by BBN
         """
         return ypBBN_to_yhe(self.Y_p(ombh2, delta_neff))
@@ -125,7 +107,7 @@ def Y_p(self, ombh2, delta_neff=0., grid=False):
         Get BBN helium nucleon fraction by intepolation in table.
 
         :param ombh2: :math:`\Omega_b h^2` (or, more generally, value of function_of[0])
-        :param delta_neff:  additional N_eff relative to standard value (of 3.046) (or value of function_of[1])
+        :param delta_neff:  additional N_eff relative to standard value (of 3.044) (or value of function_of[1])
         :param grid: parameter for :class:`~scipy:scipy.interpolate.RectBivariateSpline` (whether to evaluate the
            results on a grid spanned by the input arrays, or at points specified by the input arrays)
         :return:  Y_p helium nucleon fraction predicted by BBN. Call Y_He() to get mass fraction instead.
@@ -137,7 +119,7 @@ def DH(self, ombh2, delta_neff=0., grid=False):
         Get deuterium ratio D/H by interpolation in table
 
         :param ombh2: :math:`\Omega_b h^2` (or, more generally, value of function_of[0])
-        :param delta_neff:  additional N_eff relative to standard value (of 3.046) (or value of function_of[1])
+        :param delta_neff:  additional N_eff relative to standard value (of 3.044) (or value of function_of[1])
         :param grid: parameter for :class:`~scipy:scipy.interpolate.RectBivariateSpline` (whether to evaluate the
            results on a grid spanned by the input arrays, or at points specified by the input arrays)
         :return: D/H
@@ -151,7 +133,7 @@ def get(self, name, ombh2, delta_neff=0., grid=False):
 
         :param name: string name of the parameter, as given in header of interpolation table
         :param ombh2: :math:`\Omega_b h^2` (or, more generally, value of function_of[0])
-        :param delta_neff:  additional N_eff relative to standard value (of 3.046) (or value of function_of[1])
+        :param delta_neff:  additional N_eff relative to standard value (of 3.044) (or value of function_of[1])
         :param grid: parameter for :class:`~scipy:scipy.interpolate.RectBivariateSpline` (whether to evaluate the
            results on a grid spanned by the input arrays, or at points specified by the input arrays)
         :return:  Interpolated value (or grid)
@@ -184,23 +166,25 @@ def Y_p(self, ombh2, delta_neff=0., tau_neutron=None):
         # Parthenope fits, as in Planck 2015 papers
 
         :param ombh2: :math:`\Omega_b h^2`
-        :param delta_neff:  additional N_eff relative to standard value (of 3.046)
+        :param delta_neff:  additional N_eff relative to standard value (of 3.046 for consistency with Planck)
         :param tau_neutron: neutron lifetime
         :return:  :math:`Y_p^{\rm BBN}` helium nucleon fraction predicted by BBN
         """
+        delta_neff = default_nnu + delta_neff - 3.046
         return (
-                       0.2311 + 0.9502 * ombh2 - 11.27 * ombh2 * ombh2
-                       + delta_neff * (0.01356 + 0.008581 * ombh2 - 0.1810 * ombh2 * ombh2)
-                       + delta_neff * delta_neff * (-0.0009795 - 0.001370 * ombh2 + 0.01746 * ombh2 * ombh2)
-               ) * pow((tau_neutron or self.taun) / 880.3, 0.728)
+                0.2311 + 0.9502 * ombh2 - 11.27 * ombh2 * ombh2
+                + delta_neff * (0.01356 + 0.008581 * ombh2 - 0.1810 * ombh2 * ombh2)
+                + delta_neff * delta_neff * (-0.0009795 - 0.001370 * ombh2 + 0.01746 * ombh2 * ombh2)
+        ) * pow((tau_neutron or self.taun) / 880.3, 0.728)
 
     def DH(self, ombh2, delta_neff, tau_neutron=None):
+        delta_neff = default_nnu + delta_neff - 3.046
         return (
-                       18.754 - 1534.4 * ombh2 + 48656. * ombh2 * ombh2 - 552670. * ombh2 ** 3
-                       + delta_neff * (2.4914 - 208.11 * ombh2 + 6760.9 * ombh2 ** 2 - 78007. * ombh2 ** 3)
-                       + delta_neff * delta_neff * (
-                               0.012907 - 1.3653 * ombh2 + 37.388 * ombh2 ** 2 - 267.78 * ombh2 ** 3)
-               ) * pow((tau_neutron or self.taun) / 880.3, 0.418) * 1e-5
+                18.754 - 1534.4 * ombh2 + 48656. * ombh2 * ombh2 - 552670. * ombh2 ** 3
+                + delta_neff * (2.4914 - 208.11 * ombh2 + 6760.9 * ombh2 ** 2 - 78007. * ombh2 ** 3)
+                + delta_neff * delta_neff * (
+                        0.012907 - 1.3653 * ombh2 + 37.388 * ombh2 ** 2 - 267.78 * ombh2 ** 3)
+        ) * pow((tau_neutron or self.taun) / 880.3, 0.418) * 1e-5
 
 
 _predictors = {}

camb/camb.py

@@ -76,7 +76,7 @@ def get_age(params):
     Get age of universe for given set of parameters
 
     :param params:  :class:`.model.CAMBparams` instance
-    :return: age of universe in gigayears
+    :return: age of universe in Julian gigayears
     """
     return CAMB_GetAge(byref(params))
 

camb/constants.py

@@ -6,30 +6,31 @@
 zeta5 = 1.0369277551433699263313
 zeta7 = 1.0083492773819228268397
 
+# Constants from particle data book 2022
 c = 2.99792458e8
-h_p = 6.62606896e-34
+h_p = 6.62607015e-34
 hbar = h_p / 2 / const_pi
 
-G = 6.6738e-11  # data book 2012, last digit +/-8
-sigma_thomson = 6.6524616e-29
-sigma_boltz = 5.6704e-8
-k_B = 1.3806504e-23
-eV = 1.60217646e-19
+G = 6.67430e-11
+sigma_thomson = 6.6524587321e-29
+sigma_boltz = 5.670374419e-8
+k_B = 1.380649e-23
+eV = 1.602176634e-19
+
+m_p = 1.67262192369e-27
+m_e = 9.1093837015e-31
 
-m_p = 1.672621637e-27  # 1.672623e-27
 m_H = 1.673575e-27  # av. H atom
-m_e = 9.10938215e-31
-mass_ratio_He_H = 3.9715
+m_He4 = 6.646479073e-27  # He4
+mass_ratio_He_H = m_He4 / m_H
 
-Gyr = 3.1556926e16
-Mpc = 3.085678e22  # seem to be different definitions of this?
+Gyr = 365.25 * 86400 * 1e9
+Mpc = 3.085677581e22
 MPc_in_sec = Mpc / c  # Mpc/c = 1.029272d14 in SI units
 
 barssc0 = k_B / m_p / c ** 2
 kappa = 8. * const_pi * G
 a_rad = 8. * const_pi ** 5 * k_B ** 4 / 15 / c ** 3 / h_p ** 3
-# 7.565914e-16 #radiation constant for u=aT^4
-
 
 compton_cT = MPc_in_sec * (8. / 3.) * (sigma_thomson / (m_e * c)) * a_rad
 # compton_cT is cT in Mpc units, (8./3.)*(sigma_T/(m_e*c))*a_R in Mpc
@@ -45,7 +46,7 @@
 line21_const = 3 * l_21cm ** 2 * c * h_p / 32 / const_pi / k_B * A10 * MPc_in_sec * 1000
 # 1000 to get in MilliKelvin
 COBE_CMBTemp = 2.7255  # (Fixsen 2009) used as default value
-default_nnu = 3.046
+default_nnu = 3.044
 
 inv_neutrino_mass_fac = zeta3 * 3. / 2 / const_pi ** 2 * 4. / 11 * \
                         ((k_B * COBE_CMBTemp / hbar / c) ** 3 * kappa / 3 / (100 * 1000 / Mpc) ** 2 / (c ** 2 / eV))

camb/model.py

@@ -628,7 +628,7 @@ def get_Y_p(self, ombh2=None, delta_neff=None):
         (or the default one, if `Y_He` has not been set).
 
         :param ombh2: :math:`\Omega_b h^2` (default: value passed to :meth:`set_cosmology`)
-        :param delta_neff:  additional :math:`N_{\rm eff}` relative to standard value (of 3.046)
+        :param delta_neff:  additional :math:`N_{\rm eff}` relative to standard value (of 3.044)
                            (default: from values passed to :meth:`set_cosmology`)
         :return:  :math:`Y_p^{\rm BBN}` helium nucleon fraction predicted by BBN.
         """
@@ -646,7 +646,7 @@ def get_DH(self, ombh2=None, delta_neff=None):
         (or the default one, if `Y_He` has not been set).
 
         :param ombh2: :math:`\Omega_b h^2` (default: value passed to :meth:`set_cosmology`)
-        :param delta_neff:  additional :math:`N_{\rm eff}` relative to standard value (of 3.046)
+        :param delta_neff:  additional :math:`N_{\rm eff}` relative to standard value (of 3.044)
                            (default: from values passed to :meth:`set_cosmology`)
         :return: BBN helium nucleon fraction D/H
         """

camb/results.py

@@ -1503,7 +1503,7 @@ def hubble_parameter(self, z):
 
     def physical_time_a1_a2(self, a1, a2):
         """
-        Get physical time between two scalar factors in Gigayears
+        Get physical time between two scalar factors in Julian Gigayears
 
         Must have called :meth:`calc_background`, :meth:`calc_background_no_thermo` or calculated transfer functions
         or power spectra.
@@ -1522,7 +1522,7 @@ def physical_time_a1_a2(self, a1, a2):
 
     def physical_time(self, z):
         """
-        Get physical time from hot big bang to redshift z in Gigayears.
+        Get physical time from hot big bang to redshift z in Julian Gigayears.
 
         :param z:  redshift
         :return: t(z)/Gigayear

camb/tests/camb_test.py

@@ -124,7 +124,7 @@ def testBackground(self):
         t1 = data.conformal_time(11.5)
         t2 = data.comoving_radial_distance(11.5)
         self.assertAlmostEqual(t2, t0 - t1, 2)
-        self.assertAlmostEqual(t1, 4200.78, 2)
+        self.assertAlmostEqual(t1, 4200.809, 2)
         chistar = data.conformal_time(0) - data.tau_maxvis
         chis = np.linspace(0, chistar, 197)
         zs = data.redshift_at_comoving_radial_distance(chis)
@@ -136,25 +136,25 @@ def testBackground(self):
 
         derived = data.get_derived_params()
         self.assertAlmostEqual(derived['age'], age, 2)
-        self.assertAlmostEqual(derived['rdrag'], 146.976, 2)
+        self.assertAlmostEqual(derived['rdrag'], 146.986, 2)
         self.assertAlmostEqual(derived['rstar'], data.sound_horizon(derived['zstar']), 2)
 
         # Test BBN consistency, base_plikHM_TT_lowTEB best fit model
         pars.set_cosmology(H0=67.31, ombh2=0.022242, omch2=0.11977, mnu=0.06, omk=0)
-        self.assertAlmostEqual(pars.YHe, 0.245347, 5)
+        self.assertAlmostEqual(pars.YHe, 0.2458, 5)
         data.calc_background(pars)
-        self.assertAlmostEqual(data.cosmomc_theta(), 0.010408566, 7)
+        self.assertAlmostEqual(data.cosmomc_theta(), 0.0104090741, 7)
         self.assertAlmostEqual(data.get_derived_params()['kd'], 0.14055, 4)
 
         pars.set_cosmology(H0=67.31, ombh2=0.022242, omch2=0.11977, mnu=0.06, omk=0,
                            bbn_predictor=bbn.BBN_table_interpolator())
-        self.assertAlmostEqual(pars.YHe, 0.2453469, 5)
+        self.assertAlmostEqual(pars.YHe, 0.2458, 5)
         self.assertAlmostEqual(pars.get_Y_p(), bbn.BBN_table_interpolator().Y_p(0.022242, 0), 5)
 
         # test massive sterile models as in Planck papers
         pars.set_cosmology(H0=68.0, ombh2=0.022305, omch2=0.11873, mnu=0.06, nnu=3.073, omk=0, meffsterile=0.013)
         self.assertAlmostEqual(pars.omnuh2, 0.00078, 5)
-        self.assertAlmostEqual(pars.YHe, 0.24573, 5)
+        self.assertAlmostEqual(pars.YHe, 0.246218, 5)
         self.assertAlmostEqual(pars.N_eff, 3.073, 4)
 
         data.calc_background(pars)
@@ -175,7 +175,7 @@ def testBackground(self):
 
         # test theta
         pars.set_cosmology(cosmomc_theta=0.0104085, ombh2=0.022271, omch2=0.11914, mnu=0.06, omk=0)
-        self.assertAlmostEqual(pars.H0, 67.5512, 2)
+        self.assertAlmostEqual(pars.H0, 67.537, 2)
         with self.assertRaises(CAMBParamRangeError):
             pars.set_cosmology(cosmomc_theta=0.0204085, ombh2=0.022271, omch2=0.11914, mnu=0.06, omk=0)
         pars = camb.set_params(cosmomc_theta=0.0104077, ombh2=0.022, omch2=0.122, w=-0.95)
@@ -185,7 +185,7 @@ def testBackground(self):
         self.assertAlmostEqual(camb.get_background(pars).get_derived_params()['thetastar'] / 100, 0.010311, 7)
         pars = camb.set_params(thetastar=0.010311, ombh2=0.022, omch2=0.122, omk=-0.05)
         self.assertAlmostEqual(camb.get_background(pars).get_derived_params()['thetastar'] / 100, 0.010311, 7)
-        self.assertAlmostEqual(pars.H0, 49.7148, places=3)
+        self.assertAlmostEqual(pars.H0, 49.70624, places=3)
 
         pars = camb.set_params(cosmomc_theta=0.0104077, ombh2=0.022, omch2=0.122, w=-0.95, wa=0,
                                dark_energy_model='ppf')
@@ -328,9 +328,9 @@ def testPowers(self):
 
         data.calc_power_spectra(pars)
         kh3, z3, pk3 = data.get_matter_power_spectrum(1e-4, 1, 20)
-        self.assertAlmostEqual(pk[-1][-3], 51.909, 2)
-        self.assertAlmostEqual(pk3[-1][-3], 57.709, 2)
-        self.assertAlmostEqual(pk2[-2][-4], 56.436, 2)
+        self.assertAlmostEqual(pk[-1][-3], 51.924, 2)
+        self.assertAlmostEqual(pk3[-1][-3], 57.723, 2)
+        self.assertAlmostEqual(pk2[-2][-4], 56.454, 2)
         camb.set_feedback_level(0)
 
         PKnonlin = camb.get_matter_power_interpolator(pars, nonlinear=True)
@@ -457,10 +457,11 @@ def testSigmaR(self):
         self.assertAlmostEqual(sigma8, results.get_sigmaR(np.array([8]), z_indices=-1)[-1], places=3)
         self.assertAlmostEqual(results.get_sigmaR(8)[-1], results.get_sigmaR(8, z_indices=-1))
         pars.set_matter_power(nonlinear=False, k_per_logint=0, kmax=2)
+
         results = camb.get_results(pars)
         P, z, k = results.get_matter_power_interpolator(nonlinear=False, hubble_units=False, k_hunit=False,
                                                         return_z_k=True, extrap_kmax=100, silent=True)
-        truth = 0.800629  # from high kmax, high accuracy boost
+        truth = 0.800679  # from high kmax, high accuracy boost
         self.assertTrue(abs(results.get_sigmaR(8)[-1] / sigma8 - 1) < 1e-3)
 
         def get_sigma(_ks, dlogk):
@@ -657,7 +658,7 @@ def testSources(self):
         results = camb.get_results(pars)
         cls = results.get_source_cls_dict()
         self.assertAlmostEqual(np.sum(cls['PxW1'][10:3000:20]), 0.00020001, places=5)
-        self.assertAlmostEqual(np.sum(cls['W1xW1'][10:3000:20]), 2.26348, places=3)
+        self.assertAlmostEqual(np.sum(cls['W1xW1'][10:3000:20]), 2.26413, places=3)
         self.assertAlmostEqual(np.sum(cls['W1xW1'][10]), 0.0001097, places=6)
 
     def testSymbolic(self):

fortran/DarkEnergyFluid.f90

@@ -203,7 +203,7 @@ subroutine TAxionEffectiveFluid_Init(this, State)
             !get (very) approximate result for sound speed parameter; arXiv:1806.10608  Eq 30 (but mu may not exactly agree with what they used)
             n = nint((1+this%w_n)/(1-this%w_n))
             !Assume radiation domination, standard neutrino model; H0 factors cancel
-            grho_rad = (kappa/c**2*4*sigma_boltz/c**3*State%CP%tcmb**4*Mpc**2*(1+3.046*7._dl/8*(4._dl/11)**(4._dl/3)))
+            grho_rad = (kappa/c**2*4*sigma_boltz/c**3*State%CP%tcmb**4*Mpc**2*(1+default_nnu*7._dl/8*(4._dl/11)**(4._dl/3)))
             xc = this%a_c**2/2/sqrt(grho_rad/3)
             F=7./8
             p=1./2

fortran/camb.f90

@@ -185,7 +185,7 @@ subroutine CAMB_GetCls(State, Cls, lmax,GC_conventions)
     end subroutine CAMB_GetCls
 
     function CAMB_GetAge(P)
-    !Return age in gigayears, returns -1 on error
+    !Return age in Julian gigayears, returns -1 on error
     type(CAMBparams), intent(in) :: P
     real(dl) CAMB_GetAge
     integer error

fortran/config.f90

@@ -3,7 +3,7 @@ module config
     use constants, only: const_twopi
     implicit none
 
-    character(LEN=*), parameter :: version = '1.3.8'
+    character(LEN=*), parameter :: version = '1.4.0'
 
     integer :: FeedbackLevel = 0 !if >0 print out useful information about the model