diff --git a/.travis.yml b/.travis.yml index 7b358b5c..e0b89031 100644 --- a/.travis.yml +++ b/.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: diff --git a/camb/__init__.py b/camb/__init__.py index f41c5e66..6e077b3c 100644 --- a/camb/__init__.py +++ b/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 diff --git a/camb/bbn.py b/camb/bbn.py index 4cd210a4..48ac5aae 100644 --- a/camb/bbn.py +++ b/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 = {} diff --git a/camb/camb.py b/camb/camb.py index cca9e03a..2b9538aa 100644 --- a/camb/camb.py +++ b/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)) diff --git a/camb/cambdll.dll b/camb/cambdll.dll index 29cabf8b..dd692b9e 100644 Binary files a/camb/cambdll.dll and b/camb/cambdll.dll differ diff --git a/camb/constants.py b/camb/constants.py index 71f42789..a45b1d74 100644 --- a/camb/constants.py +++ b/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)) diff --git a/camb/model.py b/camb/model.py index b45a6b66..f47f10bf 100644 --- a/camb/model.py +++ b/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 """ diff --git a/camb/results.py b/camb/results.py index 7774cbc9..d2c5bf78 100644 --- a/camb/results.py +++ b/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 diff --git a/camb/tests/camb_test.py b/camb/tests/camb_test.py index 74bbb69f..52181633 100644 --- a/camb/tests/camb_test.py +++ b/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): diff --git a/fortran/DarkEnergyFluid.f90 b/fortran/DarkEnergyFluid.f90 index 0849e790..d6d50581 100644 --- a/fortran/DarkEnergyFluid.f90 +++ b/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 diff --git a/fortran/camb.f90 b/fortran/camb.f90 index 70acdbb8..97d08a2f 100644 --- a/fortran/camb.f90 +++ b/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 diff --git a/fortran/config.f90 b/fortran/config.f90 index bbf2d1f2..b704928e 100644 --- a/fortran/config.f90 +++ b/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 diff --git a/fortran/constants.f90 b/fortran/constants.f90 index 63a482d8..ae3dd862 100644 --- a/fortran/constants.f90 +++ b/fortran/constants.f90 @@ -19,27 +19,29 @@ module constants real(dl), parameter :: zeta5 = 1.0369277551433699263313_dl real(dl), parameter :: zeta7 = 1.0083492773819228268397_dl + !Constants from particle data book 2022 real(dl), parameter :: c = 2.99792458e8_dl - real(dl), parameter :: h_P = 6.62606896e-34_dl + real(dl), parameter :: h_P = 6.62607015e-34_dl real(dl), parameter :: hbar = h_P/const_twopi - real(dl), parameter :: G=6.6738e-11_dl !data book 2012, last digit +/-8 - real(dl), parameter :: sigma_thomson = 6.6524616e-29_dl - real(dl), parameter :: sigma_boltz = 5.6704e-8_dl - real(dl), parameter :: k_B = 1.3806504e-23_dl - real(dl), parameter :: eV = 1.60217646e-19_dl + real(dl), parameter :: G = 6.67430e-11_dl ! last two digit +/-15 + real(dl), parameter :: sigma_thomson = 6.6524587321e-29_dl + real(dl), parameter :: sigma_boltz = 5.670374419e-8_dl + real(dl), parameter :: k_B = 1.380649e-23_dl + real(dl), parameter :: eV = 1.602176634e-19_dl + real(dl), parameter :: m_p = 1.67262192369e-27_dl + real(dl), parameter :: m_e = 9.1093837015e-31_dl - real(dl), parameter :: m_p = 1.672621637e-27_dl ! 1.672623e-27_dl - real(dl), parameter :: m_H = 1.673575e-27_dl !av. H atom - real(dl), parameter :: m_e = 9.10938215e-31_dl - real(dl), parameter :: mass_ratio_He_H = 3.9715_dl + !Isotopic masses from https://www.ciaaw.org/ + real(dl), parameter :: m_H = 1.673575e-27_dl !av. H atom for D/H=2.527e-5 + real(dl), parameter :: m_He4 = 6.646479073e-27_dl + real(dl), parameter :: mass_ratio_He_H = m_He4/m_H - - real(dl), parameter :: Gyr=3.1556926e16_dl - real(dl), parameter :: Mpc = 3.085678e22_dl !seem to be different definitions of this? - real(dl), parameter :: MPC_in_sec = Mpc/c ! Mpc/c = 1.029272d14 in SI units + real(dl), parameter :: Gyr = 365.25e9_dl*86400 !Julian gigayear + real(dl), parameter :: Mpc = 3.085677581e22_dl + real(dl), parameter :: MPC_in_sec = Mpc/c ! Mpc/c = 1.029271250e14 in SI units real(dl), parameter :: barssc0= k_B / m_p / c**2 real(dl), parameter :: kappa=8._dl*const_pi*G @@ -65,7 +67,8 @@ module constants !converts omnuh2 into sum m_nu in eV for non-relativistic but thermal neutrinos (no 0.046 factor); ~ 94.07 real(dl), parameter :: neutrino_mass_fac= 1/inv_neutrino_mass_fac - real(dl), parameter :: default_nnu = 3.046_dl + !Effective number of neutrinos, arXiv:2008.01074, arXiv:2012.02726 + real(dl), parameter :: default_nnu = 3.044_dl !Neutrino mass splittings real(dl), parameter :: delta_mnu21 = 7.54e-5_dl !eV^2 Particle Data Group 2015 (-0.22, + 0.26) real(dl), parameter :: delta_mnu31 = 2.46e-3_dl !eV^2 Particle Data Group 2015 (+- 0.06) diff --git a/fortran/model.f90 b/fortran/model.f90 index 3c14cd56..7c212969 100644 --- a/fortran/model.f90 +++ b/fortran/model.f90 @@ -318,7 +318,7 @@ subroutine CAMBparams_SetNeutrinoHierarchy(this, omnuh2, omnuh2_sterile, nnu, ne neff_massive_standard=0 end if if (omnuh2_sterile>0) then - if (nnu