Repository: PlasmaPy/PlasmaPy
Pull request #957 (2020-12-04 04:27:38)

Patch mass_density() & update Alfven_speed()
This PR focuses on the `mass_density()` and `Alfven_speed()` functions.  `Alfven_speed()` is touched here because it's currently the only functionality that utilizes `mass_density()`

**`mass_density()`**

The original mass density function was a little un-intuitive since it would actually calculate the total mass density (ions + electrons) from a number density opposed to just the ion mass density or electron mass density.  The new mass density function only calculates the mass density for the particle specified.

How it now works...

```python
>>> # if a mass density is passed in then it is returned
>>> mass_density(1 * u.kg * u.m ** -3, "He")
<Quantity 1. kg/m3>

>>> # number density is assumed to be associated with the particle of interest
>>> mass_density(1e18 * u.m ** -3, "He")
<Quantity 6.64647699e-09 kg / m3>

>>> # z_ratio keyword can be used if the number density is associated with a different species
>>> ne = 1e18 * u.m ** -3
>>> par = Particle("He 2+")
>>> mass_density(ne, par, z_ratio=1/par.integer_charge)
<Quantity 3.32232756e-09 kg / m3>
```

**`Alfven_speed()`**

Original functionality does not change.  I just updated the code to handle the new `mass_density()` function and implemented handling of `Particle()`.  This `Particle` handling could be dropped/modified once the `particle_like` annotation is implemented (see issue #912).  I also restructured the tests, since that was marked as a TODO.

rocco8773: @namurphy Thanks for the review!  My code is always littered with spelling mistakes.  I swear I actively look for them before requesting a review. 🤦‍♂️🤦‍♂️🤣🤣
rocco8773: @StanczakDominik Thanks for pointing out the typos.  For some reason my IDE does not even mark them.

I also fix a few more typos I noticed while fixing your suggestions.

Patch:
diff --git a/changelog/957.breaking.rst b/changelog/957.breaking.rst
new file mode 100644
index 00000000..2583e19d
--- /dev/null
+++ b/changelog/957.breaking.rst
@@ -0,0 +1,3 @@
+Update :func:`plasmapy.formulary.parameters.mass_density` so it calculates the mass
+density for a specific particle from a given number density.  Original function
+calculated the total mass density (ion + electron).
diff --git a/changelog/957.bugfix.rst b/changelog/957.bugfix.rst
new file mode 100644
index 00000000..5704e79a
--- /dev/null
+++ b/changelog/957.bugfix.rst
@@ -0,0 +1,4 @@
+Update :func:`plasmapy.formulary.parameters.Alfven_speed` to properly use the updated
+:func:`~plasmapy.formulary.parameters.mass_density` and maintain the same behavior.
+Also add handling of the ``ion`` input keyword, so `~plasmapy.particles.Particle` and
+the `~plasmapy.particles.Particle` convertible representations can be used as inputs.
diff --git a/plasmapy/formulary/parameters.py b/plasmapy/formulary/parameters.py
index b517b720..c8f09da5 100644
--- a/plasmapy/formulary/parameters.py
+++ b/plasmapy/formulary/parameters.py
@@ -48,11 +48,12 @@
 import numpy as np
 import warnings
 
-from astropy.constants.si import c, e, eps0, k_B, m_e, m_p, mu0
-from typing import Optional
+from astropy.constants.si import c, e, eps0, k_B, mu0
+from typing import Optional, Union
 
 from plasmapy import particles
 from plasmapy.particles import Particle
+from plasmapy.particles.exceptions import ChargeError
 from plasmapy.utils import PhysicsError
 from plasmapy.utils.decorators import (
     angular_freq_to_hz,
@@ -96,64 +97,101 @@ def _grab_charge(ion: Particle, z_mean=None):
 )
 def mass_density(
     density: [u.m ** -3, u.kg / (u.m ** 3)],
-    particle: Optional[Particle] = None,
-    z_mean: Optional[numbers.Real] = None,
+    particle: Union[Particle, str],
+    z_ratio: Optional[numbers.Real] = 1,
 ) -> u.kg / u.m ** 3:
-    """Utility function to merge two possible inputs for particle charge.
+    """
+    Calculates the mass density from a number density.
+
+    .. math::
+
+        \\rho = \\left| \\frac{Z_{s}}{Z_{particle}} \\right| n_{s} m_{particle}
+              = | Z_{ratio} | n_{s} m_{particle}
+
+    where :math:`m_{particle}` is the particle mass, :math:`n_{s}` is a number
+    density for plasma species :math:`s`, :math:`Z_{s}` is the integer charge of
+    species :math:`s`, and :math:`Z_{particle}` is the integer charge of
+    ``particle``.  For example, if the electron density is given for :math:`n_s`
+    and ``particle`` is a doubly ionized atom, then :math:`Z_{ratio} = -1 / 2`.
 
     **Aliases:** `rho_`
 
     Parameters
     ----------
-    density : ~astropy.units.Quantity
-        Either a particle density (number of particles per unit volume, in units
-        of 1/m^3) or a mass density (in units of kg/m^3 or equivalent).
-
-    particle : ~plasmapy.particles.Particle, optional
-        Representation of the particle species (e.g., `'p'` for protons, `'D+'`
-        for deuterium, or `'He-4 +1'` for singly ionized helium-4).  If no
-        charge state information is provided, then the particles are assumed
-        to be singly charged. A particle is required if using a number
-        density (1/m^3) for the density parameter.
-
-    z_mean : float
-        An optional float describing the average ionization of a particle
-        species.
+    density : `~astropy.units.Quantity`
+        Either a particle number density (in units of :math:`m^{-3}` or
+        equivalent) or a mass density (in units of :math:`kg/m^3` or
+        equivalent).  If ``density`` is a mass density, then it will be passed
+        through and returned without modification.
+
+    particle : `~plasmapy.particles.Particle`
+        The particle for which the mass density is being calculated for.  Must
+        be a `~plasmapy.particles.Particle` or a value convertible to
+        a `~plasmapy.particles.Particle` (e.g., ``'p'`` for protons,
+        ``'D+'`` for deuterium, or ``'He-4 +1'`` for singly ionized helium-4).
+
+    z_ratio : `int`, `float`, optional
+        The ratio of the integer charges corresponding to the plasma species
+        represented by ``density`` and the ``particle``.  For example, if the
+        given ``density`` is and electron density and ``particle`` is doubly
+        ionized ``He``, then ``z_ratio = -0.5``.  Default is ``1``.
 
     Raises
     ------
-    ValueError
-        If the `density` has units inconvertible to either a particle density
-        or a mass density, or if you pass in a number density without a particle.
+    `~astropy.units.UnitTypeError`
+        If the ``density`` does not have units equivalent to a number density
+        or mass density.
+
+    `TypeError`
+        If ``density`` is not of type `~astropy.units.Quantity`, or convertible.
+
+    `TypeError`
+        If ``particle`` is not of type or convertible to
+        `~plasmapy.particles.Particle`.
+
+    `TypeError`
+        If ``z_ratio`` is not of type `int` or `float`.
+
+    `ValueError`
+        If ``density`` is negative.
 
     Returns
     -------
-    ~astropy.units.Quantity
-        The mass density calculated from all the provided sources of information.
+    `~astropy.units.Quantity`
+        The mass density for the plasma species represented by ``particle``.
 
     Examples
     -------
-    >>> from astropy import units as u
-    >>> mass_density(1 * u.m ** -3,'p')
-    <Quantity 1.67353...e-27 kg / m3>
-    >>> mass_density(4 * u.m ** -3,'D+')
-    <Quantity 1.33779...e-26 kg / m3>
-
+    >>> import astropy.units as u
+    >>> mass_density(1 * u.m ** -3, 'p')
+    <Quantity 1.67262...e-27 kg / m3>
+    >>> mass_density(4 * u.m ** -3, 'D+')
+    <Quantity 1.33743...e-26 kg / m3>
+    >>> mass_density(2.e12 * u.cm ** -3, 'He')
+    <Quantity 1.32929...e-08 kg / m3>
+    >>> mass_density(2.e12 * u.cm ** -3, 'He', z_ratio=0.5)
+    <Quantity 6.64647...e-09 kg / m3>
+    >>> mass_density(1.0 * u.g * u.m ** -3, "")
+    <Quantity 0.001 kg / m3>
     """
-    # validate_quantities ensures we have units of u.kg/u.m**3 or 1/u.m**3
-    rho = density
-    if not rho.unit.is_equivalent(u.kg / u.m ** 3):
-        if particle:
-            m_i = particles.particle_mass(particle)
-            Z = _grab_charge(particle, z_mean)
-            rho = density * m_i + Z * density * m_e
-        else:
-            raise ValueError(
-                f"If passing a number density, you must pass a "
-                f"particle (not {particle}) to calculate the mass density!"
+    if density.unit.is_equivalent(u.kg / u.m ** 3):
+        return density
+
+    if not isinstance(particle, Particle):
+        try:
+            particle = Particle(particle)
+        except TypeError:
+            raise TypeError(
+                f"If passing a number density, you must pass a plasmapy Particle "
+                f"(not type {type(particle)}) to calculate the mass density!"
             )
 
-    return rho
+    if not isinstance(z_ratio, (float, np.floating, int, np.integer)):
+        raise TypeError(
+            f"Expected type int or float for keyword z_ratio, got type {type(z_ratio)}."
+        )
+
+    return abs(z_ratio) * density * particle.mass
 
 
 rho_ = mass_density
@@ -166,77 +204,91 @@ def Alfven_speed(
     B: u.T,
     density: [u.m ** -3, u.kg / u.m ** 3],
     ion: Optional[Particle] = None,
-    z_mean=None,
+    z_mean: Optional[numbers.Real] = None,
 ) -> u.m / u.s:
     r"""
-    Return the Alfvén speed.
+    Calculate the Alfvén speed.
 
-    **Aliases:** `va_`
+    The Alfven speed :math:`V_A` is the typical propagation speed of magnetic
+    disturbances in a plasma, and is given by:
 
-    Parameters
-    ----------
-    B : ~astropy.units.Quantity
-        The magnetic field magnitude in units convertible to tesla.
+    .. math::
 
-    density : ~astropy.units.Quantity
-        Either the ion number density in units convertible to 1 / m**3,
-        or the mass density in units convertible to kg / m**3.
+        V_A = \frac{B}{\sqrt{\mu_0\rho}}
 
-    ion : ~plasmapy.particles.Particle, optional
-        Representation of the ion species (e.g., `'p'` for protons,
-        `'D+'` for deuterium, or `'He-4 +1'` for singly ionized
-        helium-4). If no charge state information is provided, then the
-        ions are assumed to be singly charged. If the density is an ion
-        number density, then this paramter is required in order to convert
-        to mass density.
+    where :math:`B` is the magnetic field and :math:`\rho = n_i m_i + n_e m_e`
+    is the total mass density (:math:`n_i` is the ion number density,
+    :math:`n_e` is the electron number density, :math:`m_i` is the ion mass,
+    and :math:`m_e` is the electron mass).
 
-    z_mean : ~astropy.units.Quantity, optional
-        The average ionization (arithmetic mean) for a plasma where the
-        a macroscopic description is valid. If this quantity is not
-        given then the atomic charge state (integer) of the ion
-        is used. This is effectively an average Alfven speed for the
-        plasma where multiple charge states are present.
+    **Aliases:** `va_`
+
+    Parameters
+    ----------
+    B : `~astropy.units.Quantity`
+        The magnetic field magnitude in units convertible to :math:`Tesla`.
+
+    density : `~astropy.units.Quantity`
+        Either the ion number density :math:`n_i` in units convertible to
+        :math:`m^{-3}` or the total mass density :math:`\rho` in units
+        convertible to :math:`kg / m^3`.
+
+    ion : `~plasmapy.particles.Particle`, optional
+        Representation of the ion species (e.g., `'p'` for protons, `'D+'` for
+        deuterium, `'He-4 +1'` for singly ionized helium-4, etc.). If no charge
+        state information is provided, then the ions are assumed to be singly
+        ionized. If the density is an ion number density, then this paramter
+        is required in order to convert to mass density.
+
+    z_mean : `~numbers.Real`, optional
+        The average ionization state (arithmetic mean) of the ``ion`` composing
+        the plasma.  This is used in calculating the mass density
+        :math:`\rho = n_i (m_i + Z_{mean} m_e)`.  ``z_mean`` is ignored if
+        ``density`` is passed as a mass density and overrides any charge state
+        info provided by ``ion``.
 
     Returns
     -------
-    V_A : ~astropy.units.Quantity with units of speed
-        The Alfvén speed of the plasma in units of meters per second.
+    V_A : `~astropy.units.Quantity`
+        The Alfvén speed in units :math:`m/s`.
 
     Raises
     ------
-    TypeError
-        The magnetic field and density arguments are not instances of
-        `~astropy.units.Quantity` and cannot be converted into those.
+    `~plasmapy.utils.exceptions.RelativityError`
+        If the Alfven velocity is greater than or equal to the speed of light.
 
-    ~astropy.units.UnitConversionError
-        If the magnetic field or density is not in appropriate units.
+    `TypeError`
+        If ``B`` and/or ``density`` are not of type `~astropy.units.Quantity`,
+        or convertible.
 
-    ~plasmapy.utils.RelativityError
-        If the Alfven velocity is greater than or equal to the speed of light
+    `TypeError`
+        If ``ion`` is not of type or convertible to `~plasmapy.particles.Particle`.
 
-    ValueError
-        If the density is negative, or the ion mass or charge state
-        cannot be found.
+    `TypeError`
+        If ``z_mean`` is not of type `int` or `float`.
+
+    `~astropy.units.UnitTypeError`
+        If the magnetic field ``B`` does not have units equivalent to
+        :math:`Tesla`.
+
+    `~astropy.units.UnitTypeError`
+        If the ``density`` does not have units equivalent to a number density
+        or mass density.
+
+    `ValueError`
+        If ``density`` is negative.
 
     Warns
     -----
-    ~plasmapy.utils.RelativityWarning
+    : `~plasmapy.utils.exceptions.RelativityWarning`
         If the Alfven velocity exceeds 5% of the speed of light
 
-    ~astropy.units.UnitsWarning
-        if units are not provided, SI units are assumed.
+    : `~astropy.units.UnitsWarning`
+        If units are not provided for the magnetic field ``B``, units of
+        :math:`Tesla` are assumed.
 
     Notes
     -----
-    The Alfven speed :math:`V_A` is the typical propagation speed
-    of magnetic disturbances in a plasma, and is given by:
-
-    .. math::
-
-        V_A = \frac{B}{\sqrt{\mu_0\rho}}
-
-    where the mass density is :math:`\rho = n_i m_i + n_e m_e`.
-
     This expression does not account for relativistic effects, and
     loses validity when the resulting speed is a significant fraction
     of the speed of light.
@@ -249,15 +301,39 @@ def Alfven_speed(
     >>> n = 5e19*u.m**-3
     >>> rho = n*(m_p+m_e)
     >>> ion = 'p'
-    >>> Alfven_speed(B, n, ion)
+    >>> Alfven_speed(B, n, ion=ion)
     <Quantity 43173.870... m / s>
-    >>> Alfven_speed(B, rho, ion)
+    >>> Alfven_speed(B, rho)
     <Quantity 43173.870... m / s>
-    >>> Alfven_speed(B, rho, ion).to(u.cm/u.us)
+    >>> Alfven_speed(B, rho).to(u.cm/u.us)
     <Quantity 4.317387 cm / us>
+    >>> Alfven_speed(B, n, ion="He +2")
+    <Quantity 21664.18... m / s>
+    >>> Alfven_speed(B, n, ion="He++")
+    <Quantity 21664.18... m / s>
+    >>> Alfven_speed(B, n, ion="He", z_mean=1.8)
+    <Quantity 21661.51... m / s>
 
     """
-    rho = mass_density(density, ion, z_mean)
+    if density.unit.is_equivalent(u.kg / u.m ** 3):
+        rho = density
+    else:
+        if not isinstance(ion, Particle):
+            try:
+                ion = Particle(ion)
+            except TypeError:
+                raise TypeError(
+                    f"If passing a number density, you must pass a plasmapy Particle "
+                    f"(not type {type(ion)}) to calculate the mass density!"
+                )
+        if z_mean is None:
+            try:
+                z_mean = abs(ion.integer_charge)
+            except ChargeError:
+                z_mean = 1
+
+        z_mean = abs(z_mean)
+        rho = mass_density(density, ion) + mass_density(density, "e", z_ratio=z_mean)
 
     V_A = np.abs(B) / np.sqrt(mu0 * rho)
     return V_A
@@ -1096,7 +1172,7 @@ def gyroradius(
         if isfinite_Ti:
             # T_i is valid, so use it to determine Vperp
             Vperp = thermal_speed(T_i, particle=particle)
-        # else: Vperp is alread valid, do nothing
+        # else: Vperp is already valid, do nothing
     elif np.isscalar(Vperp.value):  # only T_i is an array
         # this means either Vperp must be nan, or T_i must be array of all nan,
         # or else we couldn't have gotten through check 1


Test patch:
diff --git a/plasmapy/formulary/tests/test_parameters.py b/plasmapy/formulary/tests/test_parameters.py
index 9182de07..1efd255f 100644
--- a/plasmapy/formulary/tests/test_parameters.py
+++ b/plasmapy/formulary/tests/test_parameters.py
@@ -4,7 +4,7 @@
 import pytest
 
 from astropy import units as u
-from astropy.constants import c, e, m_e, m_p, mu0
+from astropy.constants import m_e, m_p, mu0
 from astropy.tests.helper import assert_quantity_allclose
 
 from plasmapy.formulary.parameters import (
@@ -47,6 +47,7 @@
     wp_,
     wuh_,
 )
+from plasmapy.particles import Particle
 from plasmapy.particles.exceptions import InvalidParticleError
 from plasmapy.utils.exceptions import (
     PhysicsError,
@@ -93,13 +94,50 @@
 class Test_mass_density:
     r"""Test the mass_density function in parameters.py."""
 
-    def test_particleless(self):
-        with pytest.raises(ValueError):
-            mass_density(1 * u.m ** -3)
-
-    def test_wrong_units(self):
-        with pytest.raises(u.UnitTypeError):
-            mass_density(1 * u.J)
+    @pytest.mark.parametrize(
+        "args, kwargs, conditional",
+        [
+            ((-1 * u.kg * u.m ** -3, "He"), {}, pytest.raises(ValueError)),
+            ((-1 * u.m ** -3, "He"), {}, pytest.raises(ValueError)),
+            (("not a Quantity", "He"), {}, pytest.raises(TypeError)),
+            ((1 * u.m ** -3,), {}, pytest.raises(TypeError)),
+            ((1 * u.J, "He"), {}, pytest.raises(u.UnitTypeError)),
+            ((1 * u.m ** -3, None), {}, pytest.raises(TypeError)),
+            (
+                (1 * u.m ** -3, "He"),
+                {"z_ratio": "not a ratio"},
+                pytest.raises(TypeError),
+            ),
+        ],
+    )
+    def test_raises(self, args, kwargs, conditional):
+        with conditional:
+            mass_density(*args, **kwargs)
+
+    @pytest.mark.parametrize(
+        "args, kwargs, expected",
+        [
+            ((1.0 * u.g * u.m ** -3, ""), {}, 1.0e-3 * u.kg * u.m ** -3),
+            ((5.0e12 * u.cm ** -3, "He"), {}, 3.32323849e-8 * u.kg * u.m ** -3),
+            (
+                (5.0e12 * u.cm ** -3, Particle("He")),
+                {},
+                3.32323849e-8 * u.kg * u.m ** -3,
+            ),
+            (
+                (5.0e12 * u.cm ** -3, "He"),
+                {"z_ratio": 0.5},
+                1.66161925e-08 * u.kg * u.m ** -3,
+            ),
+            (
+                (5.0e12 * u.cm ** -3, "He"),
+                {"z_ratio": -0.5},
+                1.66161925e-08 * u.kg * u.m ** -3,
+            ),
+        ],
+    )
+    def test_values(self, args, kwargs, expected):
+        assert np.isclose(mass_density(*args, **kwargs), expected)
 
     def test_handle_nparrays(self):
         """Test for ability to handle numpy array quantities"""
@@ -110,103 +148,201 @@ def test_handle_nparrays(self):
 # are generally from the NRL Plasma Formulary.
 
 
-def test_Alfven_speed():
-    r"""Test the Alfven_speed function in parameters.py."""
-
-    # TODO: break this test up until multiple tests
-
-    assert np.isclose(
-        Alfven_speed(1 * u.T, 1e-8 * u.kg * u.m ** -3, "p").value,
-        8920620.580763856,
-        rtol=1e-6,
+class TestAlfvenSpeed:
+    """Test `~plasmapy.formulary.parameters.Alfven_speed`."""
+
+    @pytest.mark.parametrize("alias", [va_])
+    def test_aliases(self, alias):
+        assert alias is Alfven_speed
+
+    @pytest.mark.parametrize(
+        "args, kwargs, _error",
+        [
+            # scenarios that raise RelativityError
+            ((10 * u.T, 1.0e-10 * u.kg * u.m ** -3), {}, RelativityError),
+            ((np.inf * u.T, 1 * u.m ** -3), {"ion": "p"}, RelativityError),
+            ((-np.inf * u.T, 1 * u.m ** -3), {"ion": "p"}, RelativityError),
+            #
+            # scenarios that raise InvalidParticleError
+            ((1 * u.T, 5e19 * u.m ** -3), {"ion": "spacecats"}, InvalidParticleError),
+            #
+            # scenarios that raise TypeError
+            (("not a Bfield", 1.0e-10 * u.kg * u.m ** -3), {}, TypeError),
+            ((10 * u.T, "not a density"), {}, TypeError),
+            ((10 * u.T, 5), {"ion": "p"}, TypeError),
+            ((1 * u.T, 1.0e18 * u.m ** -3), {"ion": ["He"]}, TypeError),
+            ((1 * u.T, 1.0e18 * u.m ** -3), {"ion": "He", "z_mean": "nope"}, TypeError),
+            #
+            # scenarios that raise UnitTypeError
+            ((1 * u.T, 1.0e18 * u.cm), {"ion": "He"}, u.UnitTypeError),
+            ((1 * u.T, 5 * u.m ** -2), {"ion": "p"}, u.UnitTypeError),
+            ((1 * u.cm, 1.0e18 * u.m ** -3), {"ion": "He"}, u.UnitTypeError),
+            ((5 * u.A, 5e19 * u.m ** -3), {"ion": "p"}, u.UnitTypeError),
+            #
+            # scenarios that raise ValueError
+            ((1 * u.T, -1.0e18 * u.m ** -3), {"ion": "He"}, ValueError),
+            (
+                (np.array([5, 6, 7]) * u.T, np.array([5, 6]) * u.m ** -3),
+                {"ion": "p"},
+                ValueError,
+            ),
+            (
+                (np.array([0.001, 0.002]) * u.T, np.array([-5e19, 6e19]) * u.m ** -3),
+                {"ion": "p"},
+                ValueError,
+            ),
+        ],
     )
+    def test_raises(self, args, kwargs, _error):
+        """Test scenarios that raise exceptions or warnings."""
+        with pytest.raises(_error):
+            Alfven_speed(*args, **kwargs)
+
+    @pytest.mark.parametrize(
+        "args, kwargs, expected, isclose_kw, _warning",
+        [
+            # scenarios that issue RelativityWarning
+            (
+                (5 * u.T, 5e19 * u.m ** -3),
+                {"ion": "H"},
+                15413707.39,
+                {},
+                RelativityWarning,
+            ),
+            (
+                (5 * u.T, 5e19 * u.m ** -3),
+                {"ion": "H+"},
+                15413707.39,
+                {"rtol": 3.0e-4},
+                RelativityWarning,
+            ),
+            (
+                (5 * u.T, 5e19 * u.m ** -3),
+                {"ion": "p"},
+                15413707.39,
+                {"rtol": 4.0e-4},
+                RelativityWarning,
+            ),
+            #
+            # scenarios that issue UnitsWarning
+            ((0.5, 1.0e18 * u.m ** -3), {"ion": "He"}, 5470657.93, {}, u.UnitsWarning),
+        ],
+    )
+    def test_warns(self, args, kwargs, expected, isclose_kw, _warning):
+        """Test scenarios that issue warnings"""
+        with pytest.warns(_warning):
+            val = Alfven_speed(*args, **kwargs)
+            assert isinstance(val, u.Quantity)
+            assert val.unit == u.m / u.s
+            assert np.isclose(val.value, expected, **isclose_kw)
+
+    @pytest.mark.parametrize(
+        "args, kwargs, expected, isclose_kw",
+        [
+            (
+                (1 * u.T, 1e-8 * u.kg * u.m ** -3),
+                {"ion": "p"},
+                8920620.58 * u.m / u.s,
+                {"rtol": 1e-6},
+            ),
+            (
+                (1 * u.T, 1e-8 * u.kg * u.m ** -3),
+                {},
+                8920620.58 * u.m / u.s,
+                {"rtol": 1e-6},
+            ),
+            (
+                (0.05 * u.T, 1e18 * u.m ** -3),
+                {"ion": "He"},
+                Alfven_speed(0.05 * u.T, 6.64738793e-09 * u.kg * u.m ** -3),
+                {},
+            ),
+            (
+                (0.05 * u.T, 1e18 * u.m ** -3),
+                {"ion": "He+"},
+                Alfven_speed(0.05 * u.T, 1e18 * u.m ** -3, ion="He"),
+                {"rtol": 7e-5},
+            ),
+            (
+                (0.05 * u.T, 1e18 * u.m ** -3),
+                {"ion": "He", "z_mean": 2},
+                Alfven_speed(0.05 * u.T, 1e18 * u.m ** -3, ion="He +2"),
+                {"rtol": 1.4e-4},
+            ),
+            (
+                (0.05 * u.T, 1e18 * u.m ** -3),
+                {"ion": Particle("He+")},
+                Alfven_speed(0.05 * u.T, 1e18 * u.m ** -3, ion="He+"),
+                {},
+            ),
+            (
+                ([0.001, 0.002] * u.T, 5e-10 * u.kg * u.m ** -3),
+                {},
+                [
+                    va_(0.001 * u.T, 5e-10 * u.kg * u.m ** -3).value,
+                    va_(0.002 * u.T, 5e-10 * u.kg * u.m ** -3).value,
+                ]
+                * (u.m / u.s),
+                {},
+            ),
+            (
+                ([0.001, 0.002] * u.T, [5e-10, 2e-10] * u.kg * u.m ** -3),
+                {},
+                [
+                    va_(0.001 * u.T, 5e-10 * u.kg * u.m ** -3).value,
+                    va_(0.002 * u.T, 2e-10 * u.kg * u.m ** -3).value,
+                ]
+                * (u.m / u.s),
+                {},
+            ),
+            (
+                (0.001 * u.T, [1.0e18, 2e18] * u.m ** -3),
+                {"ion": "p"},
+                [
+                    va_(0.001 * u.T, 1e18 * u.m ** -3, ion="p").value,
+                    va_(0.001 * u.T, 2e18 * u.m ** -3, ion="p").value,
+                ]
+                * (u.m / u.s),
+                {},
+            ),
+        ],
+    )
+    def test_values(self, args, kwargs, expected, isclose_kw):
+        """Test expected values."""
+        assert np.allclose(Alfven_speed(*args, **kwargs), expected, **isclose_kw)
+
+    @pytest.mark.parametrize(
+        "args, kwargs, nan_mask",
+        [
+            ((np.nan * u.T, 1 * u.kg * u.m ** -3), {}, []),
+            ((0.001 * u.T, np.nan * u.kg * u.m ** -3), {}, []),
+            (([np.nan, 0.001] * u.T, 1 * u.kg * u.m ** -3), {}, [True, False]),
+            (
+                (0.001 * u.T, [np.nan, 1.0, np.nan] * u.kg * u.m ** -3),
+                {},
+                [True, False, True],
+            ),
+            (([np.nan, 0.001] * u.T, [1, np.nan] * u.kg * u.m ** -3), {}, [True, True]),
+            (
+                (0.001 * u.T, [np.nan, 1e18, np.nan] * u.m ** -3),
+                {"ion": "Ar+"},
+                [True, False, True],
+            ),
+        ],
+    )
+    def test_nan_values(self, args, kwargs, nan_mask):
+        """Input scenarios that leat to `numpy.nan` values being returned."""
+        val = Alfven_speed(*args, **kwargs)
+        if np.isscalar(val.value):
+            assert np.isnan(val)
+        else:
+            nan_arr = np.isnan(val)
+            assert np.all(nan_arr[nan_mask])
+            assert np.all(not nan_arr[np.logical_not(nan_mask)])
 
-    V_A = Alfven_speed(B, n_i, "p")
-    assert np.isclose(V_A.value, (B / np.sqrt(mu0 * n_i * (m_p + m_e))).si.value)
-
-    assert Alfven_speed(B, rho, "p") == Alfven_speed(B, n_i, "p")
-
-    assert Alfven_speed(B, rho, "p").unit.is_equivalent(u.m / u.s)
-
-    assert Alfven_speed(B, rho, "p") == Alfven_speed(-B, rho, "p")
-
-    assert Alfven_speed(B, 4 * rho, "p") == 0.5 * Alfven_speed(B, rho, "p")
-
-    assert Alfven_speed(2 * B, rho, "p") == 2 * Alfven_speed(B, rho, "p")
-
-    # Case when Z=1 is assumed
-    with pytest.warns(RelativityWarning):
-        assert np.isclose(
-            Alfven_speed(5 * u.T, 5e19 * u.m ** -3, ion="H+"),
-            Alfven_speed(5 * u.T, 5e19 * u.m ** -3, ion="p"),
-            atol=0 * u.m / u.s,
-            rtol=1e-3,
-        )
-
-    # Case where magnetic field and density are Quantity arrays
-    V_A_arr = Alfven_speed(B_arr, rho_arr, "p")
-    V_A_arr0 = Alfven_speed(B_arr[0], rho_arr[0], "p")
-    V_A_arr1 = Alfven_speed(B_arr[1], rho_arr[1], "p")
-    assert np.isclose(V_A_arr0.value, V_A_arr[0].value)
-    assert np.isclose(V_A_arr1.value, V_A_arr[1].value)
-
-    # Case where magnetic field is an array but density is a scalar Quantity
-    V_A_arr = Alfven_speed(B_arr, rho, "p")
-    V_A_arr0 = Alfven_speed(B_arr[0], rho, "p")
-    V_A_arr1 = Alfven_speed(B_arr[1], rho, "p")
-    assert np.isclose(V_A_arr0.value, V_A_arr[0].value)
-    assert np.isclose(V_A_arr1.value, V_A_arr[1].value)
-
-    with pytest.raises(ValueError):
-        Alfven_speed(np.array([5, 6, 7]) * u.T, np.array([5, 6]) * u.m ** -3, "p")
-
-    assert np.isnan(Alfven_speed(B_nanarr, rho_arr, "p")[1])
-
-    with pytest.raises(ValueError):
-        Alfven_speed(B_arr, rho_negarr, "p")
-
-    with pytest.raises(u.UnitTypeError):
-        Alfven_speed(5 * u.A, n_i, ion="p")
-
-    with pytest.raises(TypeError):
-        Alfven_speed(B, 5, ion="p")
-
-    with pytest.raises(u.UnitsError):
-        Alfven_speed(B, 5 * u.m ** -2, ion="p")
-
-    with pytest.raises(InvalidParticleError):
-        Alfven_speed(B, n_i, ion="spacecats")
-
-    with pytest.warns(RelativityWarning):  # relativistic
-        Alfven_speed(5e1 * u.T, 5e19 * u.m ** -3, ion="p")
-
-    with pytest.raises(RelativityError):  # super-relativistic
-        Alfven_speed(5e8 * u.T, 5e19 * u.m ** -3, ion="p")
-
-    with pytest.raises(ValueError):
-        Alfven_speed(0.001 * u.T, -5e19 * u.m ** -3, ion="p")
-
-    assert np.isnan(Alfven_speed(np.nan * u.T, 1 * u.m ** -3, ion="p"))
-
-    assert np.isnan(Alfven_speed(1 * u.T, np.nan * u.m ** -3, ion="p"))
-
-    with pytest.raises(RelativityError):
-        assert Alfven_speed(np.inf * u.T, 1 * u.m ** -3, ion="p") == np.inf * u.m / u.s
-
-    with pytest.raises(RelativityError):
-        assert Alfven_speed(-np.inf * u.T, 1 * u.m ** -3, ion="p") == np.inf * u.m / u.s
-
-    with pytest.warns(u.UnitsWarning):
-        assert Alfven_speed(1.0, n_i, "p") == Alfven_speed(1.0 * u.T, n_i, "p")
-
-    Alfven_speed(1 * u.T, 5e19 * u.m ** -3, ion="p")
-    # testing for user input z_mean
-    testMeth1 = Alfven_speed(1 * u.T, 5e19 * u.m ** -3, ion="p", z_mean=0.8).si.value
-    testTrue1 = 3084015.75214846
-    errStr = f"Alfven_speed() gave {testMeth1}, should be {testTrue1}."
-    assert np.isclose(testMeth1, testTrue1, atol=0.0, rtol=1e-6), errStr
-
-    assert_can_handle_nparray(Alfven_speed)
+    def test_handle_nparrays(self):
+        """Test for ability to handle numpy array quantities"""
+        assert_can_handle_nparray(Alfven_speed)
 
 
 def test_ion_sound_speed():