Source code for pyccl.neutrinos

__all__ = ("NeutrinoMassSplits", "nu_masses",)

from enum import Enum
from numbers import Real
from typing import Iterable

import numpy as np
from scipy.optimize import root

from . import physical_constants as const


[docs]class NeutrinoMassSplits(Enum): """Enumeration listing all allowed neutrino mass split types. - 'sum': sum of masses. - 'single': single massive neutrino. - 'equal': total mass distributed equally among 3 species. - 'normal': normal hierarchy. - 'inverted': inverted hierarchy. - 'list': a list of 3 different masses is passed. """ SUM = 'sum' SINGLE = 'single' EQUAL = 'equal' NORMAL = 'normal' INVERTED = 'inverted' LIST = 'list' # placeholder for backwards-compatibility
[docs]def nu_masses(*, Omega_nu_h2=None, mass_split, m_nu=None): """Returns the neutrinos mass(es) for a given value of :math:`\\Omega_\\nu h^2`, according to the splitting convention specified by the user. Args: Omega_nu_h2 (:obj:`float`): Neutrino energy density at z=0 times :math:`h^2`. mass_split (:obj:`str`): indicates how the masses should be split up Should be one of 'normal', 'inverted', 'equal' or 'sum'. m_nu (:obj:`float` or array_like): Mass in eV of the massive neutrinos present. If a sequence is passed, it is assumed that the elements of the sequence represent the individual neutrino masses. Returns: :obj:`float` or `array`: Neutrino mass(es) corresponding to this :math:`\\Omega_\\nu h^2`. """ if m_nu is None: m_nu = 93.14 * Omega_nu_h2 return _get_neutrino_masses(m_nu=m_nu, mass_split=mass_split)
def _get_neutrino_masses(*, m_nu, mass_split): """ """ if isinstance(m_nu, Real) and m_nu == 0: # no massive neutrinos return np.array([]) if isinstance(m_nu, Iterable): # input was list return np.asarray(m_nu).copy() split = NeutrinoMassSplits if split(mass_split) == split.SUM: return m_nu if split(mass_split) == split.SINGLE: return np.atleast_1d(m_nu) if split(mass_split) == split.EQUAL: return np.full(3, m_nu/3) c = const D12, D13p, D13n = c.DELTAM12_sq, c.DELTAM13_sq_pos, c.DELTAM13_sq_neg def M_nu(m, D13): m2 = m * m return np.array([m.sum()-m_nu, m2[1]-m2[0]-D12, m2[2]-m2[0]-D13]) def check_mnu(val): if m_nu < val: raise ValueError(f"m_nu < {val} incompatible with mass hierarchy") if split(mass_split) == split.NORMAL: check_mnu(np.sqrt(D12) + np.sqrt(D13p)) x0 = [0, np.sqrt(D12), np.sqrt(D13p)] return root(M_nu, x0, args=(D13p,)).x if split(mass_split) == split.INVERTED: check_mnu(np.sqrt(-(D13n + D12)) + np.sqrt(-D13n)) x0 = [0, np.sqrt(-(D13n + D12)), np.sqrt(-D13n)] return root(M_nu, x0, args=(D13n,)).x