Repository: PlasmaPy/PlasmaPy

Issue:
Return namedtuples in dielectric.py functions
There are a few functions in [`plasmapy/physics/dielectric.py`](https://github.com/PlasmaPy/PlasmaPy/blob/1511167bfb7631b0b54d1e4919d0adaf001c59aa/plasmapy/physics/dielectric.py) that return three values, e.g.
```Python
def cold_plasma_permittivity_SDP(B, species, n, omega):
    ...
    return S, D, P
```
It would probably be helpful if we could return these as `collections.namedtuple` objects (or perhaps a new class) since that would group the values together in a convenient object.  If it is returned as a `namedtuple`, then it should preserve existing behavior when doing: 
```Python
>>> S, D, P = cold_plasma_permittivity_SDP(B, species, n, omega)
```
Thanks!
-Nick



Discussion:

Patch:
diff --git a/docs/development/code_guide.rst b/docs/development/code_guide.rst
index 2c84bb65..f4dc9e73 100644
--- a/docs/development/code_guide.rst
+++ b/docs/development/code_guide.rst
@@ -49,6 +49,9 @@ Coding Style
 * For multiline imports, instead of using ``\`` to break lines, wrap the
   imported names in ``()`` parentheses and use consistent whitespace.
 
+* Use ``Optional[type]`` for type hinted keyword arguments with a default value of
+  ``None``.
+
 * There should be at most one pun per 1284 lines of code.
 
 Branches, commits, and pull requests
diff --git a/plasmapy/atomic/atomic.py b/plasmapy/atomic/atomic.py
index aa9d6daa..47c54f19 100644
--- a/plasmapy/atomic/atomic.py
+++ b/plasmapy/atomic/atomic.py
@@ -1,6 +1,5 @@
 """Functions that retrieve or are related to elemental or isotopic data."""
 
-import warnings
 from typing import (
     Union,
     Optional,
@@ -11,17 +10,16 @@
 import astropy.constants as const
 import astropy.units as u
 
-from .elements import _Elements
-from .isotopes import _Isotopes
-from .particle_class import Particle
-from .particle_input import particle_input
-from .symbols import atomic_symbol
+from plasmapy.atomic.elements import _Elements
+from plasmapy.atomic.isotopes import _Isotopes
+from plasmapy.atomic.particle_class import Particle
+from plasmapy.atomic.particle_input import particle_input
+from plasmapy.atomic.symbols import atomic_symbol
 
-from ..utils import (
+from plasmapy.utils import (
     InvalidParticleError,
     InvalidElementError,
     InvalidIsotopeError,
-    InvalidIonError,
     MissingAtomicDataError,
 )
 
@@ -262,7 +260,7 @@ def particle_mass(
 
 
 @particle_input
-def isotopic_abundance(isotope: Particle, mass_numb: numbers.Integral = None) -> numbers.Real:
+def isotopic_abundance(isotope: Particle, mass_numb: Optional[numbers.Integral] = None) -> numbers.Real:
     """
     Return the isotopic abundances if known, and otherwise zero.
 
@@ -424,7 +422,7 @@ def electric_charge(particle: Particle) -> u.Quantity:
 
 
 @particle_input
-def is_stable(particle: Particle, mass_numb: numbers.Integral = None) -> bool:
+def is_stable(particle: Particle, mass_numb: Optional[numbers.Integral] = None) -> bool:
     """
     Return `True` for stable isotopes and particles and `False` for
     unstable isotopes.
@@ -478,7 +476,7 @@ def is_stable(particle: Particle, mass_numb: numbers.Integral = None) -> bool:
 
 
 @particle_input(any_of={'stable', 'unstable', 'isotope'})
-def half_life(particle: Particle, mass_numb: numbers.Integral = None) -> u.Quantity:
+def half_life(particle: Particle, mass_numb: Optional[numbers.Integral] = None) -> u.Quantity:
     """
     Return the half-life in seconds for unstable isotopes and particles,
     and numpy.inf in seconds for stable isotopes and particles.
diff --git a/plasmapy/atomic/nuclear.py b/plasmapy/atomic/nuclear.py
index 37dc75e4..2ccc0266 100644
--- a/plasmapy/atomic/nuclear.py
+++ b/plasmapy/atomic/nuclear.py
@@ -1,7 +1,7 @@
 """Functions that are related to nuclear reactions."""
 
-from astropy import units as u, constants
-from typing import List, Union
+from astropy import units as u
+from typing import List, Union, Optional
 import re
 
 from plasmapy.utils import (
@@ -15,11 +15,12 @@
 __all__ = [
     "nuclear_binding_energy",
     "nuclear_reaction_energy",
+    "mass_energy",
 ]
 
 
 @particle_input(any_of={'isotope', 'baryon'})
-def nuclear_binding_energy(particle: Particle, mass_numb: int = None) -> u.Quantity:
+def nuclear_binding_energy(particle: Particle, mass_numb: Optional[int] = None) -> u.Quantity:
     """
     Return the nuclear binding energy associated with an isotope.
 
@@ -37,7 +38,7 @@ def nuclear_binding_energy(particle: Particle, mass_numb: int = None) -> u.Quant
     Returns
     -------
     binding_energy: `~astropy.units.Quantity`
-        The binding energy of the nucleus in units of Joules.
+        The binding energy of the nucleus in units of joules.
 
     Raises
     ------
@@ -52,9 +53,12 @@ def nuclear_binding_energy(particle: Particle, mass_numb: int = None) -> u.Quant
 
     See Also
     --------
-    `~plasmapy.atomic.nuclear_reaction_energy` : Returns the change in
+    `~plasmapy.atomic.nuclear_reaction_energy` : Return the change in
         binding energy during nuclear fusion or fission reactions.
 
+    `~plasmapy.atomic.mass_energy` : Return the mass energy of a
+        nucleon or particle.
+
     Examples
     --------
     >>> from astropy import units as u
@@ -74,6 +78,46 @@ def nuclear_binding_energy(particle: Particle, mass_numb: int = None) -> u.Quant
     return particle.binding_energy.to(u.J)
 
 
+@particle_input
+def mass_energy(particle: Particle, mass_numb: Optional[int] = None) -> u.Quantity:
+    """
+    Return a particle's mass energy.  If the particle is an isotope or
+    nuclide, return the nuclear mass energy only.
+
+    Parameters
+    ----------
+    particle: `str`, `int`, or `~plasmapy.atomic.Particle`
+        A Particle object, a string representing an element or isotope,
+        or an integer representing the atomic number of an element.
+
+    mass_numb: `int` (optional)
+        The mass number of an isotope, which is required if and only
+        if the first argument can only be used to determine the
+        element and not the isotope.
+
+    Returns
+    -------
+    mass_energy: `~astropy.units.Quantity`
+        The mass energy of the particle (or in the case of ) in units of joules.
+
+    Raises
+    ------
+    `~plasmapy.utils.InvalidParticleError`
+        If the inputs do not correspond to a valid particle.
+
+    `~plasmapy.utils.AtomicError`
+        If the inputs do not correspond to a valid isotope or nucleon.
+
+    `TypeError`
+        If the inputs are not of the correct types.
+
+    >>> mass_energy('He-4')
+    <Quantity 5.97191997e-10 J>
+
+    """
+    return particle.mass_energy
+
+
 def nuclear_reaction_energy(*args, **kwargs):
     """
     Return the released energy from a nuclear reaction.
@@ -129,16 +173,19 @@ def nuclear_reaction_energy(*args, **kwargs):
     Examples
     --------
     >>> from astropy import units as u
+
     >>> nuclear_reaction_energy("D + T --> alpha + n")
     <Quantity 2.81812097e-12 J>
+
     >>> triple_alpha1 = '2*He-4 --> Be-8'
     >>> triple_alpha2 = 'Be-8 + alpha --> carbon-12'
     >>> energy_triplealpha1 = nuclear_reaction_energy(triple_alpha1)
     >>> energy_triplealpha2 = nuclear_reaction_energy(triple_alpha2)
     >>> print(energy_triplealpha1, energy_triplealpha2)
-    -1.4714307834388437e-14 J 1.18025735267267e-12 J
+    -1.4714307834388437e-14 J 1.1802573526724632e-12 J
     >>> energy_triplealpha2.to(u.MeV)
     <Quantity 7.36658704 MeV>
+
     >>> nuclear_reaction_energy(reactants=['n'], products=['p+', 'e-'])
     <Quantity 1.25343511e-13 J>
 
@@ -229,13 +276,10 @@ def add_mass_energy(particles: List[Particle]) -> u.Quantity:
         Find the total mass energy from a list of particles, while
         taking the masses of the fully ionized isotopes.
         """
-        total_mass = 0.0 * u.kg
+        total_mass_energy = 0.0 * u.J
         for particle in particles:
-            if particle.isotope:
-                total_mass += particle.nuclide_mass
-            else:
-                total_mass += particle.mass
-        return (total_mass * constants.c ** 2).to(u.J)
+            total_mass_energy += particle.mass_energy
+        return total_mass_energy.to(u.J)
 
     input_err_msg = (
         "The inputs to nuclear_reaction_energy should be either "
@@ -266,8 +310,8 @@ def add_mass_energy(particles: List[Particle]) -> u.Quantity:
 
         try:
             LHS_string, RHS_string = re.split('-+>', reaction)
-            LHS_list = re.split(' \+ ', LHS_string)
-            RHS_list = re.split(' \+ ', RHS_string)
+            LHS_list = re.split(r' \+ ', LHS_string)
+            RHS_list = re.split(r' \+ ', RHS_string)
             reactants = process_particles_list(LHS_list)
             products = process_particles_list(RHS_list)
         except Exception as ex:
diff --git a/plasmapy/atomic/parsing.py b/plasmapy/atomic/parsing.py
index c9744e3b..9ab4bee5 100644
--- a/plasmapy/atomic/parsing.py
+++ b/plasmapy/atomic/parsing.py
@@ -6,7 +6,11 @@
 from typing import (Union, Dict)
 import numbers
 
-from plasmapy.atomic.elements import (_atomic_numbers_to_symbols, _element_names_to_symbols, _Elements)
+from plasmapy.atomic.elements import (
+    _atomic_numbers_to_symbols,
+    _element_names_to_symbols,
+    _Elements)
+
 from plasmapy.atomic.isotopes import _Isotopes
 from plasmapy.atomic.special_particles import _Particles, ParticleZoo
 
diff --git a/plasmapy/atomic/particle_class.py b/plasmapy/atomic/particle_class.py
index ebf3853b..06c36379 100644
--- a/plasmapy/atomic/particle_class.py
+++ b/plasmapy/atomic/particle_class.py
@@ -276,7 +276,7 @@ class Particle:
     ``'post-transition metal'``, ``'proton'``, ``'stable'``,
     ``'transition metal'``, ``'uncharged'``, and ``'unstable'``.
 
-"""
+    """
 
     def __init__(
             self,
@@ -425,7 +425,9 @@ def __init__(
             categories.add('uncharged')
 
         if attributes['half-life'] is not None:
-            if attributes['half-life'] == np.inf * u.s:
+            if isinstance(attributes['half-life'], str):
+                categories.add('unstable')
+            elif attributes['half-life'] == np.inf * u.s:
                 categories.add('stable')
             else:
                 categories.add('unstable')
@@ -906,6 +908,129 @@ def mass(self) -> u.Quantity:
 
         raise MissingAtomicDataError(f"The mass of {self} is not available.")
 
+    @property
+    def nuclide_mass(self) -> u.Quantity:
+        """
+        Return the mass of the bare nucleus of an isotope or a neutron.
+
+        This attribute will raise a
+        `~plasmapy.utils.InvalidIsotopeError` if the particle is not an
+        isotope or neutron, or a
+        `~plasmapy.utils.MissingAtomicDataError` if the isotope mass is
+        not available.
+
+        Examples
+        --------
+        >>> deuterium = Particle('D')
+        >>> deuterium.nuclide_mass
+        <Quantity 3.34358372e-27 kg>
+
+        """
+
+        if self.isotope == 'H-1':
+            return const.m_p
+        elif self.isotope == 'D':
+            return _special_ion_masses['D 1+']
+        elif self.isotope == 'T':
+            return _special_ion_masses['T 1+']
+        elif self.particle == 'n':
+            return const.m_n
+
+        if not self.isotope:
+            raise InvalidIsotopeError(_category_errmsg(self, 'isotope'))
+
+        base_mass = self._attributes['isotope mass']
+
+        if base_mass is None:  # coverage: ignore
+            raise MissingAtomicDataError(f"The mass of a {self.isotope} nuclide is not available.")
+
+        _nuclide_mass = self._attributes['isotope mass'] - self.atomic_number * const.m_e
+
+        return _nuclide_mass.to(u.kg)
+
+    @property
+    def mass_energy(self) -> u.Quantity:
+        """
+        Return the mass energy of the particle in joules.
+
+        If the particle is an isotope or nuclide, return the mass energy
+        of the nucleus only.
+
+        If the mass of the particle is not known, then raise a
+        `~plasmapy.utils.MissingAtomicDataError`.
+
+        Examples
+        --------
+        >>> proton = Particle('p+')
+        >>> proton.mass_energy
+        <Quantity 1.50327759e-10 J>
+
+        >>> protium = Particle('H-1 0+')
+        >>> protium.mass_energy
+        <Quantity 1.50327759e-10 J>
+
+        >>> electron = Particle('electron')
+        >>> electron.mass_energy.to('MeV')
+        <Quantity 0.51099895 MeV>
+
+        """
+        try:
+            mass = self.nuclide_mass if self.isotope else self.mass
+            energy = mass * const.c ** 2
+            return energy.to(u.J)
+        except MissingAtomicDataError:
+            raise MissingAtomicDataError(
+                f"The mass energy of {self.particle} is not available "
+                f"because the mass is unknown.") from None
+
+    @property
+    def binding_energy(self) -> u.Quantity:
+        """
+        Return the nuclear binding energy in joules.
+
+        This attribute will raise an
+        `~plasmapy.utils.InvalidIsotopeError` if the particle is not a
+        nucleon or isotope.
+
+        Examples
+        --------
+        >>> alpha = Particle('alpha')
+        >>> alpha.binding_energy
+        <Quantity 4.53346938e-12 J>
+        >>> Particle('T').binding_energy.to('MeV')
+        <Quantity 8.48179621 MeV>
+
+        The binding energy of a nucleon equals 0 joules.
+
+        >>> neutron = Particle('n')
+        >>> proton = Particle('p+')
+        >>> neutron.binding_energy
+        <Quantity 0. J>
+        >>> proton.binding_energy
+        <Quantity 0. J>
+
+        """
+
+        if self._attributes['baryon number'] == 1:
+            return 0 * u.J
+
+        if not self.isotope:
+            raise InvalidIsotopeError(
+                f"The nuclear binding energy may only be calculated for nucleons and isotopes.")
+
+        number_of_protons = self.atomic_number
+        number_of_neutrons = self.mass_number - self.atomic_number
+
+        mass_of_protons = number_of_protons * const.m_p
+        mass_of_neutrons = number_of_neutrons * const.m_n
+
+        mass_of_nucleons = mass_of_protons + mass_of_neutrons
+
+        mass_defect = mass_of_nucleons - self.nuclide_mass
+        nuclear_binding_energy = mass_defect * const.c ** 2
+
+        return nuclear_binding_energy.to(u.J)
+
     @property
     def atomic_number(self) -> numbers.Integral:
         """
@@ -1086,45 +1211,6 @@ def lepton_number(self) -> numbers.Integral:
                 f"The lepton number for {self.particle} is not available.")
         return self._attributes['lepton number']
 
-    @property
-    def binding_energy(self) -> u.Quantity:
-        """
-        Return the nuclear binding energy in joules.
-
-        This attribute will raise an
-        `~plasmapy.utils.InvalidIsotopeError` if the particle is not a
-        nucleon or isotope.
-
-        Examples
-        --------
-        >>> alpha = Particle('alpha')
-        >>> alpha.binding_energy
-        <Quantity 4.53346938e-12 J>
-        >>> Particle('T').binding_energy.to('MeV')
-        <Quantity 8.48179621 MeV>
-
-        """
-
-        if self._attributes['baryon number'] == 1:
-            return 0 * u.J
-
-        if not self.isotope:
-            raise InvalidIsotopeError(
-                f"The nuclear binding energy may only be calculated for nucleons and isotopes.")
-
-        number_of_protons = self.atomic_number
-        number_of_neutrons = self.mass_number - self.atomic_number
-
-        mass_of_protons = number_of_protons * const.m_p
-        mass_of_neutrons = number_of_neutrons * const.m_n
-
-        mass_of_nucleons = mass_of_protons + mass_of_neutrons
-
-        mass_defect = mass_of_nucleons - self.nuclide_mass
-        nuclear_binding_energy = mass_defect * const.c ** 2
-
-        return nuclear_binding_energy.to(u.J)
-
     @property
     def half_life(self) -> Union[u.Quantity, str]:
         """
diff --git a/plasmapy/atomic/special_particles.py b/plasmapy/atomic/special_particles.py
index 69ab8a15..24539621 100644
--- a/plasmapy/atomic/special_particles.py
+++ b/plasmapy/atomic/special_particles.py
@@ -3,13 +3,12 @@
 information for special particles.
 """
 
-from typing import Set, Dict, List, Optional, Union
+from typing import Set, Dict
 from astropy import units as u, constants as const
 import numpy as np
 from plasmapy.atomic.elements import _PeriodicTable
 
 
-
 class _ParticleZooClass:
     """
     Create an object with taxonomy information for special particles.
diff --git a/plasmapy/atomic/symbols.py b/plasmapy/atomic/symbols.py
index 71b1fff1..25df84d3 100644
--- a/plasmapy/atomic/symbols.py
+++ b/plasmapy/atomic/symbols.py
@@ -3,6 +3,8 @@
 and numbers.
 """
 
+from typing import Optional
+
 from .particle_class import Particle
 from .particle_input import particle_input
 import numbers
@@ -94,7 +96,7 @@ def atomic_symbol(element: Particle) -> str:
 
 
 @particle_input
-def isotope_symbol(isotope: Particle, mass_numb: numbers.Integral = None) -> str:
+def isotope_symbol(isotope: Particle, mass_numb: Optional[numbers.Integral] = None) -> str:
     """
     Return the symbol representing an isotope.
 
diff --git a/plasmapy/classes/species.py b/plasmapy/classes/species.py
index a5b53842..d0686f68 100644
--- a/plasmapy/classes/species.py
+++ b/plasmapy/classes/species.py
@@ -4,7 +4,7 @@
 
 import numpy as np
 import scipy.interpolate as interp
-from ..atomic import atomic
+from plasmapy.atomic import atomic
 from astropy import constants
 from astropy import units as u
 
diff --git a/plasmapy/physics/dielectric.py b/plasmapy/physics/dielectric.py
index bf6f7e48..65042305 100644
--- a/plasmapy/physics/dielectric.py
+++ b/plasmapy/physics/dielectric.py
@@ -5,8 +5,8 @@
 from plasmapy import utils
 from plasmapy.physics import parameters
 from plasmapy.constants import (pi, m_e, c, mu0, e, eps0)
-import scipy.special as spec
 from plasmapy.mathematics import plasma_dispersion_func_deriv
+from collections import namedtuple
 
 r"""
 Values should be returned as a `~astropy.units.Quantity` in SI units.
@@ -16,6 +16,8 @@
            'cold_plasma_permittivity_LRP',
            'permittivity_1D_Maxwellian']
 
+StixTensorElements = namedtuple("StixTensorElements", ["sum", "difference", "plasma"], )
+RotatingTensorElements = namedtuple("RotatingTensorElements", ["left", "right", "plasma"], )
 
 @utils.check_quantity(
     B={'units': u.T, 'can_be_negative': False},
@@ -47,13 +49,13 @@ def cold_plasma_permittivity_SDP(B, species, n, omega):
 
     Returns
     -------
-    S : ~astropy.units.Quantity
+    sum : ~astropy.units.Quantity
         S ("Sum") dielectric tensor element.
 
-    D : ~astropy.units.Quantity
+    difference : ~astropy.units.Quantity
         D ("Difference") dielectric tensor element.
 
-    P : ~astropy.units.Quantity
+    plasma : ~astropy.units.Quantity
         P ("Plasma") dielectric tensor element.
 
     Notes
@@ -94,9 +96,11 @@ def cold_plasma_permittivity_SDP(B, species, n, omega):
     >>> species = ['e', 'D+']
     >>> n = [1e18*u.m**-3, 1e18*u.m**-3]
     >>> omega = 3.7e9*(2*pi)*(u.rad/u.s)
-    >>> S, D, P = cold_plasma_permittivity_SDP(B, species, n, omega)
+    >>> permittivity = S, D, P = cold_plasma_permittivity_SDP(B, species, n, omega)
     >>> S
     <Quantity 1.02422902>
+    >>> permittivity.sum   # namedtuple-style access
+    <Quantity 1.02422902>
     >>> D
     <Quantity 0.39089352>
     >>> P
@@ -111,7 +115,7 @@ def cold_plasma_permittivity_SDP(B, species, n, omega):
         S += - omega_p ** 2 / (omega ** 2 - omega_c ** 2)
         D += omega_c / omega * omega_p ** 2 / (omega ** 2 - omega_c ** 2)
         P += - omega_p ** 2 / omega ** 2
-    return S, D, P
+    return StixTensorElements(S, D, P)
 
 
 @utils.check_quantity(
@@ -146,13 +150,13 @@ def cold_plasma_permittivity_LRP(B: u.T, species, n, omega: u.rad / u.s):
 
     Returns
     -------
-    L : ~astropy.units.Quantity
+    left : ~astropy.units.Quantity
         L ("Left") Left-handed circularly polarization tensor element.
 
-    R : ~astropy.units.Quantity
+    right : ~astropy.units.Quantity
         R ("Right") Right-handed circularly polarization tensor element.
 
-    P : ~astropy.units.Quantity
+    plasma : ~astropy.units.Quantity
         P ("Plasma") dielectric tensor element.
 
     Notes
@@ -187,9 +191,11 @@ def cold_plasma_permittivity_LRP(B: u.T, species, n, omega: u.rad / u.s):
     >>> species = ['e', 'D+']
     >>> n = [1e18*u.m**-3, 1e18*u.m**-3]
     >>> omega = 3.7e9*(2*pi)*(u.rad/u.s)
-    >>> L, R, P = cold_plasma_permittivity_LRP(B, species, n, omega)
+    >>> L, R, P = permittivity = cold_plasma_permittivity_LRP(B, species, n, omega)
     >>> L
     <Quantity 0.63333549>
+    >>> permittivity.left    # namedtuple-style access
+    <Quantity 0.63333549>
     >>> R
     <Quantity 1.41512254>
     >>> P
@@ -204,7 +210,7 @@ def cold_plasma_permittivity_LRP(B: u.T, species, n, omega: u.rad / u.s):
         L += - omega_p ** 2 / (omega * (omega - omega_c))
         R += - omega_p ** 2 / (omega * (omega + omega_c))
         P += - omega_p ** 2 / omega ** 2
-    return L, R, P
+    return RotatingTensorElements(L, R, P)
 
 
 @u.quantity_input(omega=u.rad / u.s,
diff --git a/plasmapy/physics/parameters.py b/plasmapy/physics/parameters.py
index c5984ef6..86f78b29 100644
--- a/plasmapy/physics/parameters.py
+++ b/plasmapy/physics/parameters.py
@@ -3,6 +3,8 @@
 plasma frequency or Debye length.
 """
 
+from typing import Optional
+
 __all__ = [
     "mass_density",
     "Alfven_speed",
@@ -64,7 +66,7 @@ def _grab_charge(ion, z_mean=None):
     return Z
 
 
-def mass_density(density, particle: str = None, z_mean: numbers.Real = None) -> u.kg / u.m ** 3:
+def mass_density(density, particle: Optional[str] = None, z_mean: Optional[numbers.Real] = None) -> u.kg / u.m ** 3:
     """Utility function to merge two possible inputs for particle charge.
 
     Parameters


Test patch:
diff --git a/plasmapy/atomic/tests/test_atomic.py b/plasmapy/atomic/tests/test_atomic.py
index 4aa338f6..e3fa18ca 100644
--- a/plasmapy/atomic/tests/test_atomic.py
+++ b/plasmapy/atomic/tests/test_atomic.py
@@ -414,6 +414,7 @@ def test_periodic_table_group(self):
         with pytest.raises(TypeError):
             periodic_table_group(('B', 'Ti', 'Ge'))
 
+
 # The tables above do not include the function to be tested in order to
 # avoid cluttering up the code.  The following block of code prepends
 # the correct function to each list containing args, kwargs, and the
diff --git a/plasmapy/atomic/tests/test_nuclear.py b/plasmapy/atomic/tests/test_nuclear.py
index 9dda788a..515e374c 100644
--- a/plasmapy/atomic/tests/test_nuclear.py
+++ b/plasmapy/atomic/tests/test_nuclear.py
@@ -1,8 +1,9 @@
-from astropy import units as u
+from astropy import units as u, constants as const
 import numpy as np
-from ..nuclear import nuclear_binding_energy, nuclear_reaction_energy
+from ..nuclear import nuclear_binding_energy, nuclear_reaction_energy, mass_energy
 from ...utils import (
     InvalidParticleError,
+    InvalidIsotopeError,
     AtomicError,
     run_test,
     run_test_equivalent_calls,
@@ -17,6 +18,11 @@
     [nuclear_binding_energy, 'He-99', {}, InvalidParticleError],
     [nuclear_binding_energy, "He", {"mass_numb": 99}, InvalidParticleError],
     [nuclear_binding_energy, 3.1415926535j, {}, TypeError],
+    [mass_energy, 'e-', {}, (const.m_e * const.c ** 2).to(u.J)],
+    [mass_energy, 'p+', {}, (const.m_p * const.c ** 2).to(u.J)],
+    [mass_energy, 'H-1', {}, (const.m_p * const.c ** 2).to(u.J)],
+    [mass_energy, 'H-1 0+', {}, (const.m_p * const.c ** 2).to(u.J)],
+    [mass_energy, 'n', {}, (const.m_n * const.c ** 2).to(u.J)],
     [nuclear_reaction_energy, (), {'reactants': ['n'], 'products': 3}, TypeError],
     [nuclear_reaction_energy, (), {'reactants': ['n'], 'products': ['He-4']}, AtomicError],
     [nuclear_reaction_energy, (), {'reactants': ['h'], 'products': ['H-1']}, AtomicError],
diff --git a/plasmapy/atomic/tests/test_particle_class.py b/plasmapy/atomic/tests/test_particle_class.py
index 45256541..59b282f1 100644
--- a/plasmapy/atomic/tests/test_particle_class.py
+++ b/plasmapy/atomic/tests/test_particle_class.py
@@ -4,7 +4,7 @@
 import inspect
 from plasmapy.utils import roman
 
-from ...constants import m_p, m_e, m_n, e
+from ...constants import m_p, m_e, m_n, e, c
 
 from ...utils import (
     AtomicWarning,
@@ -140,6 +140,7 @@
       'roman_symbol': None,
       'is_ion': False,
       'mass': m_e,
+      'mass_energy': (m_e * c ** 2).to('J'),
       'nuclide_mass': InvalidIsotopeError,
       'integer_charge': 1,
       'spin': 1 / 2,
@@ -293,6 +294,7 @@
       'isotope_name': 'helium-4',
       'ionic_symbol': 'He-4 2+',
       'roman_symbol': 'He-4 III',
+      'mass_energy': 5.971919969131517e-10 * u.J,
       'is_ion': True,
       'integer_charge': 2,
       'atomic_number': 2,
@@ -303,6 +305,13 @@
       'recombine()': Particle('He-4 1+')
       }),
 
+    ('He-4 0+', {},
+     {'particle': 'He-4 0+',
+      'element': 'He',
+      'isotope': 'He-4',
+      'mass_energy': 5.971919969131517e-10 * u.J,
+      }),
+
     ('Li', {'mass_numb': 7},
      {'particle': 'Li-7',
       'element': 'Li',
@@ -335,7 +344,8 @@
       'mass_number': 276,
       'neutron_number': 164,
       'baryon_number': 276,
-      'lepton_number': 0}),
+      'lepton_number': 0,
+      }),
 
     ('muon', {},
      {'particle': 'mu-',
@@ -358,6 +368,7 @@
       'isotope': None,
       'isotope_name': InvalidElementError,
       'mass': MissingAtomicDataError,
+      'mass_energy': MissingAtomicDataError,
       'integer_charge': 0,
       'mass_number': InvalidIsotopeError,
       'element_name': InvalidElementError,
@@ -383,6 +394,8 @@
 
     (Particle('C'), {},
      {'particle': 'C',
+      'atomic_number': 6,
+      'element': 'C',
       }),
 
     (Particle('C'), {'Z': 3, 'mass_numb': 14},
diff --git a/plasmapy/atomic/tests/test_special_particles.py b/plasmapy/atomic/tests/test_special_particles.py
index 6588e28c..c2d8ca07 100644
--- a/plasmapy/atomic/tests/test_special_particles.py
+++ b/plasmapy/atomic/tests/test_special_particles.py
@@ -1,6 +1,6 @@
 import pytest
 
-from ..special_particles import (ParticleZoo, _Particles)
+from plasmapy.atomic.special_particles import (ParticleZoo, _Particles)
 
 particle_antiparticle_pairs = [
     ('e-', 'e+'),
diff --git a/plasmapy/classes/sources/tests/test_openpmd_hdf5.py b/plasmapy/classes/sources/tests/test_openpmd_hdf5.py
index 3fb4a5a4..3ce37b2c 100644
--- a/plasmapy/classes/sources/tests/test_openpmd_hdf5.py
+++ b/plasmapy/classes/sources/tests/test_openpmd_hdf5.py
@@ -16,7 +16,7 @@ class TestOpenPMD2D:
     h5 = openpmd_hdf5.HDF5Reader(hdf5=os.path.join(rootdir, "data00000255.h5"))
 
     def test_has_electric_field_with_units(self):
-        assert self.h5.electric_field.to(u.V / u.m)
+        assert isinstance(self.h5.electric_field.to(u.V / u.m), u.Quantity)
 
     def test_correct_shape_electric_field(self):
         assert self.h5.electric_field.shape == (3, 51, 201)
@@ -44,13 +44,13 @@ class TestOpenPMD3D:
     h5 = openpmd_hdf5.HDF5Reader(hdf5=os.path.join(rootdir, "data00000100.h5"))
 
     def test_has_electric_field_with_units(self):
-        assert self.h5.electric_field.to(u.V / u.m)
+        assert isinstance(self.h5.electric_field.to(u.V / u.m), u.Quantity)
 
     def test_correct_shape_electric_field(self):
         assert self.h5.electric_field.shape == (3, 26, 26, 201)
 
     def test_has_charge_density_with_units(self):
-        assert self.h5.charge_density.to(u.C / u.m**3)
+        assert isinstance(self.h5.charge_density.to(u.C / u.m**3), u.Quantity)
 
     def test_correct_shape_charge_density(self):
         assert self.h5.charge_density.shape == (26, 26, 201)
@@ -72,25 +72,25 @@ class TestOpenPMDThetaMode:
     h5 = openpmd_hdf5.HDF5Reader(hdf5=os.path.join(rootdir, "data00000200.h5"))
 
     def test_has_electric_field_with_units(self):
-        assert self.h5.electric_field.to(u.V / u.m)
+        assert isinstance(self.h5.electric_field.to(u.V / u.m), u.Quantity)
 
     def test_correct_shape_electric_field(self):
         assert self.h5.electric_field.shape == (3, 3, 51, 201)
 
     def test_has_charge_density_with_units(self):
-        assert self.h5.charge_density.to(u.C / u.m**3)
+        assert isinstance(self.h5.charge_density.to(u.C / u.m**3), u.Quantity)
 
     def test_correct_shape_charge_density(self):
         assert self.h5.charge_density.shape == (3, 51, 201)
 
     def test_has_magnetic_field_with_units(self):
-        assert self.h5.magnetic_field.to(u.T)
+        assert isinstance(self.h5.magnetic_field.to(u.T), u.Quantity)
 
     def test_correct_shape_magnetic_field(self):
         assert self.h5.magnetic_field.shape == (3, 3, 51, 201)
 
     def test_has_electric_current_with_units(self):
-        assert self.h5.electric_current.to(u.A * u.kg / u.m**3)
+        assert isinstance(self.h5.electric_current.to(u.A * u.kg / u.m**3), u.Quantity)
 
     def test_correct_shape_electric_current(self):
         assert self.h5.electric_current.shape == (3, 3, 51, 201)
diff --git a/plasmapy/physics/tests/test_dielectric.py b/plasmapy/physics/tests/test_dielectric.py
index ffcc4e16..2fb4cd9c 100644
--- a/plasmapy/physics/tests/test_dielectric.py
+++ b/plasmapy/physics/tests/test_dielectric.py
@@ -6,7 +6,9 @@
 
 from ..dielectric import (cold_plasma_permittivity_LRP,
                           cold_plasma_permittivity_SDP,
-                          permittivity_1D_Maxwellian)
+                          permittivity_1D_Maxwellian,
+                          RotatingTensorElements,
+                          StixTensorElements)
 
 from ..parameters import (plasma_frequency,
                           gyrofrequency,
@@ -46,12 +48,23 @@ def test_proton_electron_plasma(self):
         P_analytical = 1 - (omega_pe ** 2 + omega_pp ** 2) / omega ** 2
 
         species = ['e', 'p']
-        S, D, P = cold_plasma_permittivity_SDP(B, species, n, omega)
+        S, D, P = tuple_result = cold_plasma_permittivity_SDP(B, species, n, omega)
+
+        assert tuple_result.sum is S
+        assert tuple_result.difference is D
+        assert tuple_result.plasma is P
+        assert isinstance(tuple_result, StixTensorElements)
 
         assert np.isclose(S, S_analytical)
         assert np.isclose(D, D_analytical)
         assert np.isclose(P, P_analytical)
 
+        L, R, P = rotating_tuple_result = cold_plasma_permittivity_LRP(B, species, n, omega)
+        assert rotating_tuple_result.left is L
+        assert rotating_tuple_result.right is R
+        assert rotating_tuple_result.plasma is P
+        assert isinstance(rotating_tuple_result, RotatingTensorElements)
+
     def test_three_species(self):
         """
         Test with three species (2 ions): D plasma with 5%H minority fraction
@@ -79,6 +92,14 @@ def test_SD_to_LR_relationships(self):
         assert np.isclose(S, (R + L) / 2)
         assert np.isclose(D, (R - L) / 2)
 
+    def test_numpy_array_workflow(self):
+        """as per @jhillairet at https://github.com/PlasmaPy/PlasmaPy/issues/539#issuecomment-425337810 """
+        ns = np.logspace(17, 19, 50)/u.m**3
+        B0 = 4*u.T
+        omega_RF = 2*np.pi*50e6*(u.rad/u.s)
+
+        S, D, P = cold_plasma_permittivity_SDP(B=B0, species=['e', 'D+'], n=[ns, ns], omega=omega_RF)
+        assert S.shape == D.shape == P.shape == (50,)
 
 class Test_permittivity_1D_Maxwellian:
     @classmethod
diff --git a/plasmapy/physics/tests/test_parameters.py b/plasmapy/physics/tests/test_parameters.py
index ca897c70..63686127 100644
--- a/plasmapy/physics/tests/test_parameters.py
+++ b/plasmapy/physics/tests/test_parameters.py
@@ -884,7 +884,7 @@ def test_magnetic_energy_density():
 
     assert_quantity_allclose(magnetic_energy_density(B), magnetic_energy_density(B.to(u.G)))
 
-    assert magnetic_energy_density(B_arr)
+    assert isinstance(magnetic_energy_density(B_arr), u.Quantity)
 
     with pytest.warns(u.UnitsWarning):
         magnetic_energy_density(5)