'Neutrino density function relative to the energy density in
photons.
Parameters
z : array like
Redshift
Returns
f : ndarray, or float if z is scalar
The neutrino density scaling factor relative to the density
in photons at each redshift
Notes
The density in neutrinos is given by
.. math::
\rho_{\nu} \left(a\right) = 0.2271 \, N_{eff} \,
f\left(m_{\nu} a / T_{\nu 0} \right) \,
\rho_{\gamma} \left( a \right)
where
.. math::
f \left(y\right) = \frac{120}{7 \pi^4}
\int_0^{\infty} \, dx \frac{x^2 \sqrt{x^2 + y^2}}
{e^x + 1}
assuming that all neutrino species have the same mass.
If they have different masses, a similar term is calculated
for each one. Note that f has the asymptotic behavior :math:`f(0) = 1`.
This method returns :math:`0.2271 f` using an
analytical fitting formula given in Komatsu et al. 2011, ApJS 192, 18.'

def nu_relative_density(self, z):

    prefac = 0.22710731766
    if (not self._massivenu):
        if np.isscalar(z):
            return (prefac * self._Neff)
        else:
            return ((prefac * self._Neff) * np.ones(np.asanyarray(z).shape, dtype=np.float))
    p = 1.83
    invp = 0.54644808743
    k = 0.3173
    z = np.asarray(z)
    curr_nu_y = (self._nu_y / (1.0 + np.expand_dims(z, axis=(-1))))
    rel_mass_per = ((1.0 + ((k * curr_nu_y) ** p)) ** invp)
    rel_mass = (rel_mass_per.sum((-1)) + self._nmasslessnu)
    return ((prefac * self._neff_per_nu) * rel_mass)
