Strange behaviour in classy for massless neutrinos

[arthurmloureiro]

Hi!

I am doing some tests in the neutrino sector of classy for using the pk_cb and I noticed some inconsistent behaviour.

Case 1: N_Ncdm = 0

In this case, if any m_ncdm is given, classy will raise an exception saying that m_ncdm was not used.

class_fiducial_params = {
    [...]
    'N_ncdm': 0,  # No massive neutrinos
}

# Linear Model 1
results_class_M1 = Class()
results_class_M1.set(class_fiducial_params)
results_class_M1.compute()

print(results_class_M1.Om_ncdm(0.0))
>>>  np.float64(0.0)

Otherwise, the 0.0 value is expected.

Case 2: N_ncdm = 1, m_ncdm =0.0

In this case, I would expect getting an exception or zero neutrino energy density, instead I get:

# Linear power spectrum parameters for CLASS
class_fiducial_params = {
    [...]
    'm_ncdm': 0.0, # Still no massive neutrinos
    'N_ncdm': 1,  
}

# Linear Model 1
results_class_M1 = Class()
results_class_M1.set(class_fiducial_params)
results_class_M1.compute()

print(results_class_M1.Om_ncdm(0.0))
>>> 1.2418798823407054e-05

Now this is different than what we see in case 1.

Is this behaviour expected or is this a bug?