Fix omeganuh2 (#808)

* remove massless neutrino catch

* fix flake8

* reincorporate massless option

pyccl/ccl_neutrinos.i

@@ -14,8 +14,8 @@
     if numpy.shape(a) != (nout,):
         raise CCLError("Input shape for `a` must match `(nout,)`!")
 
-    if numpy.shape(mnu) != (3,):
-        raise CCLError("Input shape for `mnu` must match `(3,)`!")
+    if numpy.shape(mnu) != (N_nu_mass,):
+        raise CCLError("Input shape for `mnu` must match `(N_nu_mass,)`!")
 %}
 
 %inline %{

pyccl/neutrinos.py

@@ -13,18 +13,16 @@
 
 def Omeganuh2(a, m_nu, T_CMB=None):
     """Calculate :math:`\\Omega_\\nu\\,h^2` at a given scale factor given
-    the sum of the neutrino masses.
-
-    .. note:: for all practical purposes, :math:`N_{\\rm eff}` is simply
-        `N_nu_mass`.
+    the neutrino masses.
 
     Args:
         a (float or array-like): Scale factor, normalized to 1 today.
-        m_nu (float or array-like): Neutrino mass (in eV)
+        m_nu (float or array-like): Neutrino mass(es) (in eV)
         T_CMB (float, optional): Temperature of the CMB (K). Default: 2.725.
 
     Returns:
-        float or array_like: corresponding to a given neutrino mass.
+        float or array_like: :math:`\\Omega_\\nu\\,h^2` at a given
+        scale factor given the neutrino masses
     """
     status = 0
     scalar = True if np.ndim(a) == 0 else False

src/ccl_neutrinos.c

@@ -170,14 +170,7 @@ double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double T_CMB, int* s
 
   Tnu = T_CMB*pow(4./11.,1./3.);
   a4 = a*a*a*a;
-  // Check if mnu=0. We assume that in the massless case mnu is a
-  // pointer to a single element and that element is 0. This should in principle never be called.
-  if (mnu[0] < 0.00017) {  // Limit taken from Lesgourges et al. 2012
-    prefix_massless = NU_CONST  * Tnu * Tnu * Tnu * Tnu;
-    return N_nu_mass*prefix_massless*7./8./a4;
-  }
-
-  // And the remaining massive case.
+
   // Tnu_eff is used in the massive case because CLASS uses an effective
   // temperature of nonLCDM components to match to mnu / Omeganu =93.14eV. Tnu_eff = T_ncdm * T_CMB = 0.71611 * T_CMB
   Tnu_eff = Tnu * ccl_constants.TNCDM / (pow(4./11.,1./3.));
@@ -187,13 +180,22 @@ double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double T_CMB, int* s
 
   OmNuh2 = 0.; // Initialize to 0 - we add to this for each massive neutrino species.
   for(int i=0; i < N_nu_mass; i++) {
-    // Get mass over T (mass (eV) / ((kb eV/s/K) Tnu_eff (K))
-    // This returns the density normalized so that we get nuh2 at a=0
-    mnuOT = mnu[i] / (Tnu_eff/a) * (ccl_constants.EV_IN_J / (ccl_constants.KBOLTZ));
-
-    // Get the value of the phase-space integral
-    intval = nu_phasespace_intg(mnuOT, status);
-    OmNuh2 = intval*prefix_massive/a4 + OmNuh2;
+
+    // Check whether this species is effectively massless
+    // In this case, invoke the analytic massless limit:
+    if (mnu[i] < 0.00017) {  // Limit taken from Lesgourges et al. 2012
+      prefix_massless = NU_CONST  * Tnu * Tnu * Tnu * Tnu;	
+      OmNuh2 = N_nu_mass*prefix_massless*7./8./a4 + OmNuh2;	
+    } else {  
+       // For the true massive case:  
+       // Get mass over T (mass (eV) / ((kb eV/s/K) Tnu_eff (K))
+       // This returns the density normalized so that we get nuh2 at a=0
+       mnuOT = mnu[i] / (Tnu_eff/a) * (ccl_constants.EV_IN_J / (ccl_constants.KBOLTZ));
+
+       // Get the value of the phase-space integral
+       intval = nu_phasespace_intg(mnuOT, status);
+       OmNuh2 = intval*prefix_massive/a4 + OmNuh2;
+    }
   }
 
   return OmNuh2;