Division by zero in BeamCXLine._beam_population(), BeamEmissionLine._beam_emission_rate and SingleRayAttenuator_beam_stopping() if plasma contains neutral atoms.

[vsnever]

Here is an example from SingleRayAttenuator: https://github.com/cherab/core/blob/2cf200ae209c11d00863602c89a5b6bb60deb85e/cherab/core/model/attenuator/singleray.pyx#L252-L268

For neutral atoms, target_z is zero and target_equiv_ne is inf. We don't get errors here because the stopping rate data for neutral atoms is missing and NullBeamStoppingRate() is used. However, I think division by zero should be avoided anyway. BeamCXLine and BeamEmissionLine have similar issue.

Since this stopping model is applicable for charge particles only I suggest adding only charged particles in the _stopping_data:

        self._stopping_data = []
        for species in self._plasma.composition:
            if species.charge:
                stopping_coeff = self._atomic_data.beam_stopping_rate(self._beam.element, species.element, species.charge)
                self._stopping_data.append((species, stopping_coeff))

It will also remove unnecessary calculations for neutrals. @mateasek, @jacklovell, what's your opinion on this?

[CnlPepper]

This calculation (derived from the adas docs) is only intended for ionised species. As Adas will never supply rates for neutrals a check should be unnecessary.

You really should not make an implicit choice for the user by just ignoring species, it would be best to raise an error that the model does not support neutral species. The user ahould then be able to configure the model to ignore neutrals. Remember cherab is meant to show all choices explicitly to avoid user error/generating results that are unexpected.

[vsnever]

An error will be raised if OpenADAS is initialised with missing_rates_return_null=False (default). However, in the simulations with multiple plasma models (e.g. active and passive CX), some plasma models may require neutral atoms, and the neutral atoms must be added to plasma compositions. In this case missing_rates_return_null must be set to True, which will cause the neutral atoms to have NullBeamStoppingRate(), and that's why the code will not produce any errors. But we will still have these useless calculations and division by zero target_z.

What I mean is that if the model works only for ionised species, neutral species should simply be ignored in this model, even if they are present in the plasma. If we want an error to be raised in _populate_stopping_data_cache() if missing_rates_return_null is False, then let's do the check in _beam_stopping():

        for species, coeff in self._stopping_data:

            if species.charge == 0:
                continue

            target_z = species.charge
            target_ne = species.distribution.density(x, y, z) * target_z
            target_ti = species.distribution.effective_temperature(x, y, z)
            target_velocity = species.distribution.bulk_velocity(x, y, z)

[vsnever]

Looks like it's a duplicate of #149, so I'm closing it. After reading the discussion in #149, I also think the neutral atoms should not be excluded until a more complex model is implemented. If the user will provide the stopping rates for neutrals, returning nan is better than giving the false impression that the user data has been accounted for in the model.