Repository: PlasmaPy/PlasmaPy

Issue:
Importing a plasmapy module fails if h5py isn't installed
```python
    from plasmapy.physics import dimensionless
  File "/Users/dstansby/github/plasmapy/plasmapy/__init__.py", line 29, in <module>
    from . import classes
  File "/Users/dstansby/github/plasmapy/plasmapy/classes/__init__.py", line 4, in <module>
    from . import sources
  File "/Users/dstansby/github/plasmapy/plasmapy/classes/sources/__init__.py", line 3, in <module>
    from .openpmd_hdf5 import HDF5Reader
  File "/Users/dstansby/github/plasmapy/plasmapy/classes/sources/openpmd_hdf5.py", line 1, in <module>
    import h5py
ModuleNotFoundError: No module named 'h5py'
```
The problem is that `__init__.py` in the `plasmapy` directory tries to import *everything*, which means e.g. `h5py` is required for plasmapy to work. I'm guessing this is an optional dependency, and shouldn't be required?

Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/561?src=pr&el=h1) Report
> Merging [#561](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/561?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/ec90f19064721cfae35edd63c2330c26d9bb46f7?src=pr&el=desc) will **decrease** coverage by `<.01%`.
> The diff coverage is `n/a`.

[![Impacted file tree graph](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/561/graphs/tree.svg?width=650&token=LzZ5mnXVad&height=150&src=pr)](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/561?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master     #561      +/-   ##
==========================================
- Coverage   97.44%   97.43%   -0.01%     
==========================================
  Files          46       46              
  Lines        3835     3827       -8     
==========================================
- Hits         3737     3729       -8     
  Misses         98       98
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/561?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/561/diff?src=pr&el=tree#diff-cGxhc21hcHkvX19pbml0X18ucHk=) | `56.25% <ø> (-14.59%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/561?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/561?src=pr&el=footer). Last update [ec90f19...9e14346](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/561?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

dstansby: I think it might be most useful to explicitly specify optional dependencies for each sub-module of plasmapy, instead of for certain functionality. That way `pip install plasmapy[submodule]` is a nice way of doing "I want to use this sub-module, this command will install all the dependencies I need".
StanczakDominik: Ooooh. That's a really neat idea!
StanczakDominik: Okay, now this seems rather silly. Apparently, some of our .rst documentation files were not getting checked for doctests (I have no idea what caused them to get tested now). The first of these is not a big deal - a nuclear reaction syntax example.

But the second is `plasmapy.test()`.

In a doctest.

And this is getting launched now.
StanczakDominik: I did forget to update the docs to account for @dstansby 's idea, so that's one thing!
StanczakDominik: All in all, I'm not sure the split of optional dependencies that we have right now (as per David's suggestion) is optimal. I have no idea how to describe this in docs. What we might need to do is reorganize our subpackages again:
* merge `math`, `physics` and `transport` into, maybe, `theory` (this would then involve lmfit and mpmath)
* keep diagnostics depending on `matplotlib` (maybe even throw in an optional pandas interface because if you're working with experimental data you most likely want something compatible with dataframes)
* rename `classes` somehow because that name is completely out  of whack right now (there's a bunch of classes in each module now). I'm not sure what a good name would be.
* remove unused `io` and `visualization` subpackages, maybe even get rid of `plasmapy.constants` because they're wrappers to astropy anyway.

But that's probably out of scope for this PR. For optional dependency name sets, I'm kind of pondering going back to what we had before 414c and leaving the above changes for a future PR.

Thoughts?

Patch:
diff --git a/.coveragerc b/.coveragerc
index 61ea1967..332aba95 100644
--- a/.coveragerc
+++ b/.coveragerc
@@ -8,6 +8,5 @@ omit =
     plasmapy/version.py
 exclude_lines =
     coverage: ignore
-    def __repr__
     ImportError
     ModuleNotFoundError
diff --git a/docs/install.rst b/docs/install.rst
index 7cd479bf..0bcd2291 100644
--- a/docs/install.rst
+++ b/docs/install.rst
@@ -22,12 +22,18 @@ installation:
 - `NumPy <http://www.numpy.org/>`_ 1.13 or newer
 - `SciPy <https://www.scipy.org/>`_ 0.19 or newer
 - `Astropy <http://www.astropy.org/>`_ 2.0 or newer
-- `matplotlib <https://matplotlib.org/>`_ 2.0 or newer
 - `Cython <http://cython.org/>`_ 0.27.2 or newer
+- `colorama <https://pypi.org/project/colorama/>`_ 0.3 or newer
+
+PlasmaPy also uses the following optional dependencies:
+
+- `matplotlib <https://matplotlib.org/>`_ 2.0 or newer
 - `h5py <https://www.h5py.org/>`_ 2.8 or newer
 - `mpmath <http://mpmath.org/>`_ 1.0 or newer
 - `lmfit <https://lmfit.github.io/lmfit-py/>`_ 0.9.7 or newer
-- `colorama <https://pypi.org/project/colorama/>`_ 0.3 or newer
+
+PlasmaPy can be installed with all of the optional dependencies via
+`pip install plasmapy[optional]`.
 
 .. _create-conda-env:
 
@@ -71,12 +77,20 @@ Installation with pip
 ---------------------
 
 To install the most recent release of PlasmaPy `on PyPI`_ with `pip
-<https://pip.pypa.io/en/stable/>`_ into an existing Python environment,
-run
+<https://pip.pypa.io/en/stable/>`_ into an existing Python environment
+alongside all optional dependencies, run
+
+.. code:: bash
+
+   pip install plasmapy[optional]
+
+To install a minimal set of dependencies (which does not guarantee that
+everything will run and may result in `ImportError`s, skip `[all]` and run
+simply
 
 .. code:: bash
 
-   pip install plasmapy
+    pip install plasmapy
 
 .. _install-conda:
 
@@ -85,6 +99,8 @@ Installation with conda
 
 We’re not on conda just yet, but we’re working on it!
 
+Conda installs all dependencies by default.
+
 Building and installing from source code
 ========================================
 
diff --git a/plasmapy/__init__.py b/plasmapy/__init__.py
index be054fa4..23eabcc3 100644
--- a/plasmapy/__init__.py
+++ b/plasmapy/__init__.py
@@ -23,16 +23,6 @@ class UnsupportedPythonError(Exception):
 if sys.version_info < tuple((int(val) for val in "3.6".split('.'))):
     raise UnsupportedPythonError("plasmapy does not support Python < {}".format(3.6))
 
-if not _ASTROPY_SETUP_:
-    # For egg_info test builds to pass, put package imports here.
-    from . import atomic
-    from . import classes
-    from . import constants
-    from . import diagnostics
-    from . import mathematics
-    from . import physics
-    from . import utils
-
 def online_help(query):
     """
     Search the online PlasmaPy documentation for the given query from plasmapy.org
diff --git a/plasmapy/atomic/__init__.py b/plasmapy/atomic/__init__.py
index f9fc6d7a..872e7810 100644
--- a/plasmapy/atomic/__init__.py
+++ b/plasmapy/atomic/__init__.py
@@ -37,6 +37,11 @@
     nuclear_reaction_energy,
 )
 
+from .ionization_state import IonizationState, State
+from .ionization_states import IonizationStates
+
+# Create instances of the most commonly used particles
+
 proton = Particle("p+")
 electron = Particle("e-")
 neutron = Particle("n")
diff --git a/plasmapy/atomic/atomic.py b/plasmapy/atomic/atomic.py
index 47c54f19..c89ce4fa 100644
--- a/plasmapy/atomic/atomic.py
+++ b/plasmapy/atomic/atomic.py
@@ -7,6 +7,7 @@
     Any,
 )
 
+import numpy as np
 import astropy.constants as const
 import astropy.units as u
 
diff --git a/plasmapy/atomic/ionization_state.py b/plasmapy/atomic/ionization_state.py
new file mode 100644
index 00000000..ae833765
--- /dev/null
+++ b/plasmapy/atomic/ionization_state.py
@@ -0,0 +1,665 @@
+"""
+Objects for storing ionization state data for a single element or for
+a single ionization level.
+"""
+
+from numbers import Integral, Real
+from typing import Union, List, Optional
+import collections
+import warnings
+
+import numpy as np
+import astropy.units as u
+
+from plasmapy.atomic import Particle, particle_input
+from plasmapy.utils import AtomicError, ChargeError, InvalidParticleError, check_quantity
+
+__all__ = ["IonizationState", "State"]
+
+
+_number_density_errmsg = (
+    "Number densities must be Quantity objects with units of inverse "
+    "volume."
+)
+
+# TODO: Change `State` into a class with validations for all of the
+# TODO: attributes.
+
+State = collections.namedtuple(
+    'State', [
+        'integer_charge',
+        'ionic_fraction',
+        'ionic_symbol',
+        'number_density',
+    ])
+
+
+class IonizationState:
+    """
+    Representation of the ionization state distribution of a single
+    element or isotope.
+
+    Parameters
+    ----------
+    particle: str, integer, or ~plasmapy.atomic.Particle
+        A `str` or `~plasmapy.atomic.Particle` instance representing
+        an element or isotope, or an integer representing the atomic
+        number of an element.
+
+    ionic_fractions: ~numpy.ndarray, list, tuple, or ~astropy.units.Quantity; optional
+        The ionization fractions of an element, where the indices
+        correspond to integer charge.  This argument should contain the
+        atomic number plus one items, and must sum to one within an
+        absolute tolerance of ``tol`` if dimensionless.  Alternatively,
+        this argument may be a `~astropy.units.Quantity` that represents
+        the number densities of each neutral/ion.
+
+    T_e: ~astropy.units.Quantity, keyword-only, optional
+        The electron temperature or thermal energy per particle.
+
+    n_elem: ~astropy.units.Quantity, keyword-only, optional
+        The number density of the element, including neutrals and all
+        ions.
+
+    tol: float or integer, keyword-only, optional
+        The absolute tolerance used by `~numpy.isclose` when testing
+        normalizations and making comparisons.  Defaults to ``1e-15``.
+
+    Raises
+    ------
+    ~plasmapy.utils.AtomicError
+        If the ionic fractions are not normalized or contain invalid
+        values, or if number density information is provided through
+        both ``ionic_fractions`` and ``n_elem``.
+
+    ~plasmapy.utils.InvalidParticleError
+        If the particle is invalid.
+
+    Examples
+    --------
+    >>> states = IonizationState('H', [0.6, 0.4], n_elem=1*u.cm**-3, T_e=11000*u.K)
+    >>> states.ionic_fractions[0]  # fraction of hydrogen that is neutral
+    0.6
+    >>> states.ionic_fractions[1]  # fraction of hydrogen that is ionized
+    0.4
+    >>> states.n_e  # electron number density
+    <Quantity 400000. 1 / m3>
+    >>> states.n_elem  # element number density
+    <Quantity 1000000. 1 / m3>
+
+    Notes
+    -----
+    Calculation of collisional ionization equilibrium has not yet been
+    implemented.
+
+    """
+
+    # TODO: Allow this class to (optionally?) handle negatively charged
+    # TODO: ions.  There are instances where singly negatively charged
+    # TODO: ions are important in astrophysical plasmas, such as H- in
+    # TODO: the atmospheres of relatively cool stars.  There may be some
+    # TODO: rare situations where doubly negatively charged ions show up
+    # TODO: too, but triply negatively charged ions are very unlikely.
+
+    # TODO: Add in functionality to find equilibrium ionization states.
+
+    @check_quantity(
+        T_e={"units": u.K},
+        n_elem={"units": u.m ** -3},
+    )
+    @particle_input(require='element', exclude='ion')
+    def __init__(self,
+                 particle: Particle,
+                 ionic_fractions=None,
+                 *,
+                 T_e: u.K = np.nan * u.K,
+                 kappa: Real = np.inf,
+                 n_elem: u.m ** -3 = np.nan * u.m ** -3,
+                 tol: Union[float, int] = 1e-15):
+        """Initialize an `~plasmapy.atomic.IonizationState` instance."""
+
+        self._particle_instance = particle
+
+        try:
+            self.tol = tol
+            self.T_e = T_e
+            self.kappa = kappa
+
+            if not np.isnan(n_elem) and isinstance(ionic_fractions, u.Quantity) and \
+                    ionic_fractions.si.unit == u.m ** -3:
+                raise AtomicError(
+                    "Cannot simultaneously provide number density "
+                    "through both n_elem and ionic_fractions.")
+
+            self.n_elem = n_elem
+            self.ionic_fractions = ionic_fractions
+
+            if ionic_fractions is None and not np.isnan(self.T_e):
+                warnings.warn(
+                    "Collisional ionization equilibration has not yet "
+                    "been implemented in IonizationState; cannot set "
+                    "ionic fractions.")
+
+        except Exception as exc:
+            raise AtomicError(
+                f"Unable to create IonizationState instance for "
+                f"{particle.particle}.") from exc
+
+    def __str__(self) -> str:
+        return f"<IonizationState instance for {self.base_particle}>"
+
+    def __repr__(self) -> str:
+        return self.__str__()
+
+    def __getitem__(self, value) -> State:
+        """Return information for a single ionization level."""
+        if isinstance(value, slice):
+            raise TypeError("IonizationState instances cannot be sliced.")
+
+        if isinstance(value, Integral) and 0 <= value <= self.atomic_number:
+            result = State(
+                value,
+                self.ionic_fractions[value],
+                self.ionic_symbols[value],
+                self.number_densities[value],
+            )
+        else:
+            if not isinstance(value, Particle):
+                try:
+                    value = Particle(value)
+                except InvalidParticleError as exc:
+                    raise InvalidParticleError(
+                        f"{value} is not a valid integer charge or "
+                        f"particle.") from exc
+
+            same_element = value.element == self.element
+            same_isotope = value.isotope == self.isotope
+            has_charge_info = value.is_category(any_of=["charged", "uncharged"])
+
+            if same_element and same_isotope and has_charge_info:
+                Z = value.integer_charge
+                result = State(
+                    Z,
+                    self.ionic_fractions[Z],
+                    self.ionic_symbols[Z],
+                    self.number_densities[Z],
+                )
+            else:
+                if not same_element or not same_isotope:
+                    raise AtomicError("Inconsistent element or isotope.")
+                elif not has_charge_info:
+                    raise ChargeError("No integer charge provided.")
+        return result
+
+    def __setitem__(self, key, value):
+        raise NotImplementedError(
+            "Item assignment of an IonizationState instance is not "
+            "allowed because the ionic fractions for different "
+            "ionization levels must be set simultaneously due to the "
+            "normalization constraint.")
+
+    def __iter__(self):
+        """Initialize an instance prior to iteration."""
+        self._charge_index = 0
+        return self
+
+    def __next__(self):
+        """
+        Return a `~plasmapy.atomic.State` instance that contains
+        information about a particular ionization level.
+        """
+        if self._charge_index <= self.atomic_number:
+            result = State(
+                self._charge_index,
+                self._ionic_fractions[self._charge_index],
+                self.ionic_symbols[self._charge_index],
+                self.number_densities[self._charge_index],
+            )
+            self._charge_index += 1
+            return result
+        else:
+            del self._charge_index
+            raise StopIteration
+
+    def __eq__(self, other):
+        """
+        Return `True` if the ionic fractions, number density scaling
+        factor (if set), and electron temperature (if set) are all
+        equal, and `False` otherwise.
+
+        Raises
+        ------
+        TypeError
+            If ``other`` is not an `~plasmapy.atomic.IonizationState`
+            instance.
+
+        AtomicError
+            If ``other`` corresponds to a different element or isotope.
+
+        Examples
+        --------
+        >>> IonizationState('H', [1, 0], tol=1e-6) == IonizationState('H', [1, 1e-6], tol=1e-6)
+        True
+        >>> IonizationState('H', [1, 0], tol=1e-8) == IonizationState('H', [1, 1e-6], tol=1e-5)
+        False
+
+        """
+        if not isinstance(other, IonizationState):
+            raise TypeError(
+                "An instance of the IonizationState class may only be "
+                "compared with another IonizationState instance.")
+
+        same_element = self.element == other.element
+        same_isotope = self.isotope == other.isotope
+
+        if not same_element or not same_isotope:
+            raise AtomicError(
+                "An instance of the IonizationState class may only be "
+                "compared with another IonizationState instance if "
+                "both correspond to the same element and/or isotope.")
+
+        # Use the tighter of the two tolerances. For thermodynamic
+        # quantities, use it as a relative tolerance because the values
+        # may substantially depart from order unity.
+
+        min_tol = np.min([self.tol, other.tol])
+
+        same_T_e = np.isnan(self.T_e) and np.isnan(other.T_e) or \
+            u.allclose(self.T_e, other.T_e, rtol=min_tol*u.K, atol=0*u.K)
+
+        same_n_elem = np.isnan(self.n_elem) and np.isnan(other.n_elem) or \
+            u.allclose(self.n_elem, other.n_elem, rtol=min_tol*u.m**-3, atol=0*u.m**-3)
+
+        # For the next line, recall that np.nan == np.nan is False (sigh)
+
+        same_fractions = np.any([
+            np.allclose(self.ionic_fractions, other.ionic_fractions, rtol=0, atol=min_tol),
+            np.all(np.isnan(self.ionic_fractions)) and np.all(np.isnan(other.ionic_fractions)),
+        ])
+
+        return np.all([same_element, same_isotope, same_T_e, same_n_elem, same_fractions])
+
+    @property
+    def ionic_fractions(self) -> np.ndarray:
+        """
+        Return the ionic fractions, where the index corresponds to
+        the integer charge.
+
+        Examples
+        --------
+        >>> hydrogen_states = IonizationState('H', [0.9, 0.1])
+        >>> hydrogen_states.ionic_fractions
+        array([0.9, 0.1])
+
+        """
+        return self._ionic_fractions
+
+    @ionic_fractions.setter
+    def ionic_fractions(self, fractions):
+        """
+        Set the ionic fractions, while checking that the new values are
+        valid and normalized to one.
+        """
+        if fractions is None or np.all(np.isnan(fractions)):
+            self._ionic_fractions = np.full(self.atomic_number + 1, np.nan, dtype=np.float64)
+            return
+
+        try:
+            if np.min(fractions) < 0:
+                raise AtomicError("Cannot have negative ionic fractions.")
+
+            if len(fractions) != self.atomic_number + 1:
+                raise AtomicError(
+                    "The length of ionic_fractions must be "
+                    f"{self.atomic_number + 1}.")
+
+            if isinstance(fractions, u.Quantity):
+                fractions = fractions.to(u.m ** -3)
+                self.n_elem = np.sum(fractions)
+                self._ionic_fractions = np.array(fractions/self.n_elem)
+            else:
+                fractions = np.array(fractions, dtype=np.float64)
+                sum_of_fractions = np.sum(fractions)
+                all_nans = np.all(np.isnan(fractions))
+
+                if not all_nans:
+                    if np.any(fractions < 0) or np.any(fractions > 1):
+                        raise AtomicError("Ionic fractions must be between 0 and 1.")
+
+                    if not np.isclose(sum_of_fractions, 1, rtol=0, atol=self.tol):
+                        raise AtomicError("Ionic fractions must sum to one.")
+
+                self._ionic_fractions = fractions
+
+        except Exception as exc:
+            raise AtomicError(f"Unable to set ionic fractions of {self.element} "
+                              f"to {fractions}.") from exc
+
+    def _is_normalized(self, tol: Optional[Real] = None) -> bool:
+        """
+        Return `True` if the sum of the ionization fractions is equal to
+        one within the allowed tolerance, and `False` otherwise.
+        """
+        tol = tol if tol is not None else self.tol
+        if not isinstance(tol, Real):
+            raise TypeError("tol must be an int or float.")
+        if not 0 <= tol < 1:
+            raise ValueError("Need 0 <= tol < 1.")
+        total = np.sum(self._ionic_fractions)
+        return np.isclose(total, 1, atol=tol, rtol=0)
+
+    def normalize(self) -> None:
+        """
+        Normalize the ionization state distribution (if set) so that the
+        sum becomes equal to one.
+        """
+        self._ionic_fractions = self._ionic_fractions / np.sum(self._ionic_fractions)
+
+    @property
+    def equil_ionic_fractions(self, T_e: u.K = None):
+        """
+        Return the equilibrium ionic fractions for temperature ``T_e``
+        or the temperature set in the IonizationState instance.  Not
+        implemented.
+        """
+        raise NotImplementedError
+
+    @u.quantity_input(equivalencies=u.temperature_energy())
+    def equilibrate(self, T_e: u.K = np.nan * u.K):
+        """
+        Set the ionic fractions to collisional ionization equilibrium
+        for temperature ``T_e``.  Not implemented.
+        """
+        # self.ionic_fractions = self.equil_ionic_fractions
+        raise NotImplementedError
+
+    @property
+    @u.quantity_input
+    def n_e(self) -> u.m ** -3:
+        """
+        Return the electron number density assuming a single species
+        plasma.
+        """
+        return np.sum(self._n_elem * self.ionic_fractions * self.integer_charges)
+
+    @property
+    @u.quantity_input
+    def n_elem(self) -> u.m ** -3:
+        """Return the total number density of neutrals and all ions."""
+        return self._n_elem.to(u.m ** -3)
+
+    @n_elem.setter
+    @u.quantity_input
+    def n_elem(self, value: u.m ** -3):
+        """Set the number density of neutrals and all ions."""
+        if value < 0 * u.m ** -3:
+            raise AtomicError
+        if 0 * u.m ** -3 < value <= np.inf * u.m ** -3:
+            self._n_elem = value.to(u.m ** -3)
+        elif np.isnan(value):
+            self._n_elem = np.nan * u.m ** -3
+
+    @property
+    @u.quantity_input
+    def number_densities(self) -> u.m ** -3:
+        """Return the number densities for each state."""
+        try:
+            return (self.n_elem * self.ionic_fractions).to(u.m ** -3)
+        except Exception:
+            return np.full(self.atomic_number + 1, np.nan) * u.m ** -3
+
+    @number_densities.setter
+    @u.quantity_input
+    def number_densities(self, value: u.m ** -3):
+        """Set the number densities for each state."""
+        if np.any(value.value < 0):
+            raise AtomicError("Number densities cannot be negative.")
+        if len(value) != self.atomic_number + 1:
+            raise AtomicError(
+                f"Incorrect number of charge states for "
+                f"{self.base_particle}")
+        value = value.to(u.m ** -3)
+
+        self._n_elem = value.sum()
+        self._ionic_fractions = value / self._n_elem
+
+    @property
+    @u.quantity_input(equivalencies=u.temperature_energy())
+    def T_e(self) -> u.K:
+        """Return the electron temperature."""
+        if self._T_e is None:
+            raise AtomicError("No electron temperature has been specified.")
+        return self._T_e.to(u.K, equivalencies=u.temperature_energy())
+
+    @T_e.setter
+    @u.quantity_input(equivalencies=u.temperature_energy())
+    def T_e(self, value: u.K):
+        """Set the electron temperature."""
+        try:
+            value = value.to(u.K, equivalencies=u.temperature_energy())
+        except (AttributeError, u.UnitsError, u.UnitConversionError):
+            raise AtomicError("Invalid temperature.") from None
+        else:
+            if value < 0 * u.K:
+                raise AtomicError("T_e cannot be negative.")
+        self._T_e = value
+
+    @property
+    def kappa(self) -> np.real:
+        """
+        Return the kappa parameter for a kappa distribution function
+        for electrons.
+
+        The value of ``kappa`` must be greater than ``1.5`` in order to
+        have a valid distribution function.  If ``kappa`` equals
+        `~numpy.inf`, then the distribution function reduces to a
+        Maxwellian.
+
+        """
+        return self._kappa
+
+    @kappa.setter
+    def kappa(self, value: Real):
+        """
+        Set the kappa parameter for a kappa distribution function for
+        electrons.  The value must be between ``1.5`` and `~numpy.inf`.
+        """
+        kappa_errmsg = "kappa must be a real number greater than 1.5"
+        if not isinstance(value, Real):
+            raise TypeError(kappa_errmsg)
+        if value <= 1.5:
+            raise ValueError(kappa_errmsg)
+        self._kappa = np.real(value)
+
+    @property
+    def element(self) -> str:
+        """Return the atomic symbol of the element."""
+        return self._particle_instance.element
+
+    @property
+    def isotope(self) -> Optional[str]:
+        """
+        Return the isotope symbol for an isotope, or `None` if the
+        particle is not an isotope.
+        """
+        return self._particle_instance.isotope
+
+    @property
+    def base_particle(self) -> str:
+        """Return the symbol of the element or isotope."""
+        return self.isotope if self.isotope else self.element
+
+    @property
+    def atomic_number(self) -> int:
+        """Return the atomic number of the element."""
+        return self._particle_instance.atomic_number
+
+    @property
+    def _particle_instances(self) -> List[Particle]:
+        """
+        Return a list of the `~plasmapy.atomic.Particle` class
+        instances corresponding to each ion.
+        """
+        return [
+            Particle(self._particle_instance.particle, Z=i)
+            for i in range(self.atomic_number + 1)
+        ]
+
+    @property
+    def ionic_symbols(self) -> List[str]:
+        """Return the ionic symbols for all charge states."""
+        return [particle.ionic_symbol for particle in self._particle_instances]
+
+    @property
+    def integer_charges(self) -> np.ndarray:
+        """Return an array with the integer charges."""
+        return np.arange(0, self.atomic_number + 1, dtype=np.int)
+
+    @property
+    def Z_mean(self) -> np.float64:
+        """Return the mean integer charge"""
+        if np.nan in self.ionic_fractions:
+            raise ChargeError(
+                "Z_mean cannot be found because no ionic fraction "
+                f"information is available for {self.base_particle}.")
+        return np.sum(self.ionic_fractions * np.arange(self.atomic_number + 1))
+
+    @property
+    def Z_rms(self) -> np.float64:
+        """Return the root mean square integer charge."""
+        return np.sqrt(np.sum(self.ionic_fractions * np.arange(self.atomic_number + 1) ** 2))
+
+    @property
+    def Z_most_abundant(self) -> List[Integral]:
+        """
+        Return a `list` of the integer charges with the highest ionic
+        fractions.
+
+        Examples
+        --------
+        >>> He = IonizationState('He', [0.2, 0.5, 0.3])
+        >>> He.Z_most_abundant
+        [1]
+        >>> Li = IonizationState('Li', [0.4, 0.4, 0.2, 0.0])
+        >>> Li.Z_most_abundant
+        [0, 1]
+
+        """
+        if np.any(np.isnan(self.ionic_fractions)):
+            raise AtomicError(
+                f"Cannot find most abundant ion of {self.base_particle} "
+                f"because the ionic fractions have not been defined.")
+
+        return np.flatnonzero(
+            self.ionic_fractions == self.ionic_fractions.max()
+        ).tolist()
+
+    @property
+    def tol(self) -> Real:
+        """Return the absolute tolerance for comparisons."""
+        return self._tol
+
+    @tol.setter
+    def tol(self, atol: Real):
+        """Set the absolute tolerance for comparisons."""
+        if not isinstance(atol, Real):
+            raise TypeError("The attribute tol must be a real number.")
+        if 0 <= atol < 1:
+            self._tol = atol
+        else:
+            raise ValueError("Need 0 <= tol < 1.")
+
+
+    def _get_states_info(self, minimum_ionic_fraction=0.01) -> List[str]:
+        """
+        Return a `list` containing the ion symbol, ionic fraction, and
+        (if available) the number density for that ion.
+        """
+
+        states_info = []
+
+        for state in self:
+            if state.ionic_fraction > minimum_ionic_fraction:
+                state_info = ""
+                symbol = state.ionic_symbol
+                if state.integer_charge < 10:
+                    symbol = symbol[:-2] + ' ' + symbol[-2:]
+                fraction = "{:.3f}".format(state.ionic_fraction)
+
+                state_info += f'{symbol}: {fraction}'
+
+                if np.isfinite(self.n_elem):
+                    value = "{:.2e}".format(state.number_density.si.value)
+                    state_info += f"    n_i = {value} m**-3"
+
+                states_info.append(state_info)
+
+        return states_info
+
+    def info(self, minimum_ionic_fraction: Real = 0.01) -> None:
+        """
+        Print quicklook information for an
+        `~plasmapy.atomic.IonizationState` instance.
+
+        Parameters
+        ----------
+        minimum_ionic_fraction: Real
+            If the ionic fraction for a particular ionization state is
+            below this level, then information for it will not be
+            printed.  Defaults to 0.01.
+
+        Example
+        -------
+        >>> He_states = IonizationState(
+        ...     'He',
+        ...     [0.941, 0.058, 0.001],
+        ...     T_e = 5.34 * u.K,
+        ...     kappa = 4.05,
+        ...     n_elem = 5.51e19 * u.m ** -3,
+        ... )
+        >>> He_states.info()
+        IonizationState instance for He with Z_mean = 0.06
+        ----------------------------------------------------------------
+        He  0+: 0.941    n_i = 5.18e+19 m**-3
+        He  1+: 0.058    n_i = 3.20e+18 m**-3
+        ----------------------------------------------------------------
+        n_elem = 5.51e+19 m**-3
+        n_e = 3.31e+18 m**-3
+        T_e = 5.34e+00 K
+        kappa = 4.05
+        ----------------------------------------------------------------
+
+        """
+        separator_line = [64 * '-']
+
+        scientific = "{:.2e}"
+        floaty = "{:.2f}"
+
+        n_elem = scientific.format(self.n_elem.value)
+        n_e = scientific.format(self.n_e.value)
+        T_e = scientific.format(self.T_e.value)
+        kappa = floaty.format(self.kappa)
+        Z_mean = floaty.format(self.Z_mean)
+
+        output = [f"IonizationState instance for {self.base_particle} with Z_mean = {Z_mean}"]
+        attributes = []
+
+        if not np.all(np.isnan(self.ionic_fractions)):
+            output += separator_line
+            output += self._get_states_info(minimum_ionic_fraction)
+            output += separator_line
+
+        if not np.isnan(self.n_elem):
+            attributes.append(f"n_elem = {n_elem} m**-3")
+            attributes.append(f"n_e = {n_e} m**-3")
+        if not np.isnan(self.T_e):
+            attributes.append(f"T_e = {T_e} K")
+        if np.isfinite(self.kappa):
+            attributes.append(f"kappa = {kappa}")
+
+        if attributes:
+            attributes += separator_line
+            output += attributes
+
+        for line in output:
+            print(line)
diff --git a/plasmapy/atomic/ionization_states.py b/plasmapy/atomic/ionization_states.py
new file mode 100644
index 00000000..a414c559
--- /dev/null
+++ b/plasmapy/atomic/ionization_states.py
@@ -0,0 +1,918 @@
+"""
+A class for storing ionization state data for multiple elements or
+isotopes.
+"""
+
+from numbers import Real, Integral
+from typing import Dict, List, Optional, Tuple, Union
+import collections
+
+import numpy as np
+import astropy.units as u
+
+from plasmapy.atomic import atomic_number, Particle, particle_symbol, IonizationState, State
+from plasmapy.utils import AtomicError, ChargeError, InvalidParticleError, check_quantity
+
+__all__ = ["IonizationStates"]
+
+
+class IonizationStates:
+    """
+    Describe the ionization state distributions of multiple elements
+    or isotopes.
+
+    Parameters
+    ----------
+    inputs: `list`, `tuple`, or `dict`
+        A `list` or `tuple` of elements or isotopes (if ``T_e`` is
+        provided); a `list` of `~plasmapy.atomic.IonizationState`
+        instances; a `dict` with elements or isotopes as keys and
+        a `~numpy.ndarray` of ionic fractions as the values; or a `dict`
+        with elements or isotopes as keys and `~astropy.units.Quantity`
+        instances with units of number density.
+
+    abundances: `dict` or `str`, optional, keyword-only
+        The relative abundances of each element in the plasma.
+
+    log_abundances: `dict`, optional, keyword-only
+        The base 10 logarithm of the relative abundances of each element
+        in the plasma.
+
+    n: ~astropy.units.Quantity, optional, keyword-only
+        The number density scaling factor.  The number density of an
+        element will be the product of its abundance and ``n``.
+
+    T_e: `~astropy.units.Quantity`, optional, keyword-only
+        The electron temperature in units of temperature or thermal
+        energy per base_particle.
+
+    kappa: float, optional, keyword-only
+        The value of kappa for a kappa distribution function.
+
+    tol: float or integer, keyword-only, optional
+        The absolute tolerance used by `~numpy.isclose` when testing
+        normalizations and making comparisons.  Defaults to ``1e-15``.
+
+    equilibrate: `bool`, optional, keyword-only
+        Set the ionic fractions to the estimated collisional ionization
+        equilibrium.  Not implemented.
+
+    Raises
+    ------
+    AtomicError
+        If `~plasmapy.atomic.IonizationStates` cannot be instantiated.
+
+    Examples
+    --------
+    >>> from astropy import units as u
+    >>> from plasmapy.atomic import IonizationStates
+    >>> states = IonizationStates(
+    ...     {'H': [0.5, 0.5], 'He': [0.95, 0.05, 0]},
+    ...     T_e = 1.2e4 * u.K,
+    ...     n = 1e15 * u.m ** -3,
+    ...     abundances = {'H': 1, 'He': 0.08},
+    ... )
+    >>> states.ionic_fractions
+    {'H': array([0.5, 0.5]), 'He': array([0.95, 0.05, 0.  ])}
+
+    The number densities are given by the ionic fractions multiplied by
+    the abundance and the
+
+    >>> states.number_densities['H']
+    <Quantity [5.e+14, 5.e+14] 1 / m3>
+    >>> states['He'] = [0.4, 0.59, 0.01]
+
+    To change the ionic fractions for a single element, use item
+    assignment.
+
+    >>> states = IonizationStates(['H', 'He'])
+    >>> states['H'] = [0.1, 0.9]
+
+    Item assignment will also work if you supply number densities.
+
+    >>> states['He'] = [0.4, 0.6, 0.0] * u.m ** -3
+    >>> states.ionic_fractions['He']
+    array([0.4, 0.6, 0. ])
+    >>> states.number_densities['He']
+    <Quantity [0.4, 0.6, 0. ] 1 / m3>
+
+    Notes
+    -----
+    No more than one of ``abundances``, ``log_abundances``, and
+    ``number_densities`` may be specified.
+
+    If the value provided during item assignment is a
+    `~astropy.units.Quantity` with units of number density that retains
+    the total element density, then the ionic fractions will be set
+    proportionately.
+
+    When making comparisons between `~plasmapy.atomic.IonizationStates`
+    instances, `~numpy.nan` values are treated as equal.  Equality tests
+    are performed to within a tolerance of ``tol``.
+
+    Collisional ionization equilibrium is based on atomic data that
+    has relative errors of order 20%.
+
+    """
+
+    # TODO: The docstring above needs to be expanded and revised to
+    # TODO: better describe what the magic methods do.
+
+    @check_quantity(
+        T_e={"units": u.K},
+        n={"units": u.m ** -3},
+    )
+    def __init__(
+            self,
+            inputs: Union[Dict[str, np.ndarray], List, Tuple],
+            *,
+            T_e: u.K = np.nan * u.K,
+            equilibrate: Optional[bool] = None,
+            abundances: Optional[Dict[str, Real]] = None,
+            log_abundances: Optional[Dict[str, Real]] = None,
+            n: u.m ** -3 = np.nan * u.m ** -3,
+            tol: Real = 1e-15,
+            kappa: Real = np.inf):
+
+        abundances_provided = abundances is not None or log_abundances is not None
+
+        set_abundances = True
+        if isinstance(inputs, dict):
+            all_quantities = np.all([isinstance(fracs, u.Quantity) for fracs in inputs.values()])
+            if all_quantities:
+                right_units = np.all([fracs[0].si.unit == u.m ** -3 for fracs in inputs.values()])
+                if not right_units:
+                    raise AtomicError("Units must be inverse volume for number densities.")
+                if abundances_provided:
+                    raise AtomicError(
+                        "Abundances cannot be provided if inputs "
+                        "provides number density information.")
+                set_abundances = False
+
+        try:
+            self._pars = collections.defaultdict(lambda: None)
+            self.T_e = T_e
+            self.n = n
+            self.tol = tol
+            self.ionic_fractions = inputs
+            if set_abundances:
+                self.abundances = abundances
+                self.log_abundances = log_abundances
+            self.kappa = kappa
+        except Exception as exc:
+            raise AtomicError("Unable to create IonizationStates instance.") from exc
+
+        if equilibrate:
+            self.equilibrate()  # for now, this raises a NotImplementedError
+
+    def __str__(self) -> str:
+        return f"<IonizationStates for: {', '.join(self.base_particles)}>"
+
+    def __repr__(self) -> str:
+        return self.__str__()
+
+    def __getitem__(self, *values) -> IonizationState:
+
+        errmsg = f"Invalid indexing for IonizationStates instance: {values[0]}"
+
+        one_input = not isinstance(values[0], tuple)
+        two_inputs = len(values[0]) == 2
+
+        if not one_input and not two_inputs:
+            raise IndexError(errmsg)
+
+        try:
+            arg1 = values[0] if one_input else values[0][0]
+            int_charge = None if one_input else values[0][1]
+            particle = arg1 if arg1 in self.base_particles else particle_symbol(arg1)
+
+            if int_charge is None:
+                return IonizationState(
+                    particle=particle,
+                    ionic_fractions=self.ionic_fractions[particle],
+                    T_e=self._pars["T_e"],
+                    n_elem=np.sum(self.number_densities[particle]),
+                    tol=self.tol,
+                )
+            else:
+                if not isinstance(int_charge, Integral):
+                    raise TypeError(f"{int_charge} is not a valid charge for {base_particle}.")
+                elif not 0 <= int_charge <= atomic_number(particle):
+                    raise ChargeError(f"{int_charge} is not a valid charge for {base_particle}.")
+                return State(
+                    integer_charge=int_charge,
+                    ionic_fraction=self.ionic_fractions[particle][int_charge],
+                    ionic_symbol=particle_symbol(particle, Z=int_charge),
+                    number_density=self.number_densities[particle][int_charge]
+                )
+        except Exception as exc:
+            raise IndexError(errmsg) from exc
+
+    def __setitem__(self, key, value):
+
+        errmsg = (
+            f"Cannot set item for this IonizationStates instance for "
+            f"key = {repr(key)} and value = {repr(value)}")
+
+        try:
+            particle = particle_symbol(key)
+            self.ionic_fractions[key]
+        except (AtomicError, TypeError):
+            raise KeyError(
+                f"{errmsg} because {repr(key)} is an invalid particle."
+            ) from None
+        except KeyError:
+            raise KeyError(
+                f"{errmsg} because {repr(key)} is not one of the base "
+                f"particles whose ionization state is being kept track "
+                f"of.") from None
+
+        if isinstance(value, u.Quantity) and value.unit != u.dimensionless_unscaled:
+            try:
+                new_number_densities = value.to(u.m ** -3)
+            except u.UnitConversionError:
+                raise ValueError(
+                    f"{errmsg} because the units of value do not "
+                    f"correspond to a number density.") from None
+
+            old_n_elem = np.sum(self.number_densities[particle])
+            new_n_elem = np.sum(new_number_densities)
+
+            density_was_nan = np.all(np.isnan(self.number_densities[particle]))
+            same_density = u.quantity.allclose(old_n_elem, new_n_elem, rtol=self.tol)
+
+            if not same_density and not density_was_nan:
+                raise ValueError(
+                    f"{errmsg} because the old element number density "
+                    f"of {old_n_elem} is not approximately equal to "
+                    f"the new element number density of {new_n_elem}.")
+
+            value = (new_number_densities / new_n_elem).to(u.dimensionless_unscaled)
+
+            # If the abundance of this particle has not been defined,
+            # then set the abundance if there is enough (but not too
+            # much) information to do so.
+
+            abundance_is_undefined = np.isnan(self.abundances[particle])
+            isnan_of_abundance_values = np.isnan(list(self.abundances.values()))
+            all_abundances_are_nan = np.all(isnan_of_abundance_values)
+            n_is_defined = not np.isnan(self.n)
+
+            if abundance_is_undefined:
+                if n_is_defined:
+                    self._pars['abundances'][particle] = new_n_elem / self.n
+                elif all_abundances_are_nan:
+                    self.n = new_n_elem
+                    self._pars['abundances'][particle] = 1
+                else:
+                    raise AtomicError(
+                        f"Cannot set number density of {particle} to "
+                        f"{value * new_n_elem} when the number density "
+                        f"scaling factor is undefined, the abundance "
+                        f"of {particle} is undefined, and some of the "
+                        f"abundances of other elements/isotopes is "
+                        f"defined.")
+
+        try:
+            new_fractions = np.array(value, dtype=np.float64)
+        except Exception as exc:
+            raise TypeError(
+                f"{errmsg} because value cannot be converted into an "
+                f"array that represents ionic fractions.") from exc
+
+        # TODO: Create a separate function that makes sure ionic
+        # TODO: fractions are valid to reduce code repetition.  This
+        # TODO: would probably best go as a private function in
+        # TODO: ionization_state.py.
+
+        required_nstates = atomic_number(particle) + 1
+        new_nstates = len(new_fractions)
+        if new_nstates != required_nstates:
+            raise ValueError(
+                f"{errmsg} because value must have {required_nstates} "
+                f"ionization levels but instead corresponds to "
+                f"{new_nstates} levels.")
+
+        all_nans = np.all(np.isnan(new_fractions))
+        if not all_nans and (new_fractions.min() < 0 or new_fractions.max() > 1):
+            raise ValueError(
+                f"{errmsg} because the new ionic fractions are not "
+                f"all between 0 and 1.")
+
+        normalized = np.isclose(np.sum(new_fractions), 1, rtol=self.tol)
+        if not normalized and not all_nans:
+            raise ValueError(
+                f"{errmsg} because the ionic fractions are not "
+                f"normalized to one.")
+
+        self._ionic_fractions[particle][:] = new_fractions[:]
+
+    def __iter__(self):
+        """
+        Prepare an `~plasmapy.atomic.IonizationStates` instance for
+        iteration.
+        """
+        self._element_index = 0
+        return self
+
+    @property
+    def __ITER__(self):  # coverage: ignore
+        """
+        Recall that our code development guide states that there should
+        be at most one pun per 1284 lines of code.
+        """
+        raise NotImplementedError(
+            "The International Thermonuclear Experimental Reactor "
+            "is still under construction.")
+
+    def __next__(self):
+        if self._element_index < len(self.base_particles):
+            particle = self.base_particles[self._element_index]
+            result = IonizationState(
+                particle,
+                self.ionic_fractions[particle],
+                T_e=self.T_e,
+                n_elem=np.sum(self.number_densities[particle]),
+                tol=self.tol,
+            )
+            self._element_index += 1
+            return result
+        else:
+            del self._element_index
+            raise StopIteration
+
+    def __eq__(self, other):
+
+        if not isinstance(other, IonizationStates):
+            raise TypeError(
+                "IonizationStates instance can only be compared with "
+                "other IonizationStates instances.")
+
+        if self.base_particles != other.base_particles:
+            raise AtomicError(
+                "Two IonizationStates instances can be compared only "
+                "if the base particles are the same.")
+
+        min_tol = np.min([self.tol, other.tol])
+
+        # Check any of a whole bunch of equality measures, recalling
+        # that np.nan == np.nan is False.
+
+        for attribute in ['T_e', 'n_e', 'kappa']:
+            this = eval(f"self.{attribute}")
+            that = eval(f"other.{attribute}")
+
+            # TODO: Maybe create a function in utils called same_enough
+            # TODO: that would take care of all of these disparate
+            # TODO: equality measures.
+
+            this_equals_that = np.any([
+                this == that,
+                this is that,
+                np.isnan(this) and np.isnan(that),
+                np.isinf(this) and np.isinf(that),
+                u.quantity.allclose(this, that, rtol=min_tol),
+            ])
+
+            if not this_equals_that:
+                return False
+
+        for attribute in ['ionic_fractions', 'number_densities']:
+
+            this_dict = eval(f"self.{attribute}")
+            that_dict = eval(f"other.{attribute}")
+
+            for particle in self.base_particles:
+
+                this = this_dict[particle]
+                that = that_dict[particle]
+
+                this_equals_that = np.any([
+                    this is that,
+                    np.all(np.isnan(this)) and np.all(np.isnan(that)),
+                    u.quantity.allclose(this, that, rtol=min_tol),
+                ])
+
+                if not this_equals_that:
+                    return False
+
+        return True
+
+    @property
+    def ionic_fractions(self) -> Dict[str, np.array]:
+        """
+        Return a `dict` containing the ionic fractions for each element
+        and isotope.
+
+        The keys of this `dict` are the symbols for each element or
+        isotope.  The values will be `~numpy.ndarray` objects containing
+        the ionic fractions for each ionization level corresponding to
+        each element or isotope.
+
+        """
+        return self._ionic_fractions
+
+    @ionic_fractions.setter
+    def ionic_fractions(self, inputs: Union[Dict, List, Tuple]):
+        """
+        Set the ionic fractions.
+
+        Notes
+        -----
+        The ionic fractions are initialized during instantiation of
+        `~plasmapy.atomic.IonizationStates`.  After this, the only way
+        to reset the ionic fractions via the ``ionic_fractions``
+        attribute is via a `dict` with elements or isotopes that are a
+        superset of the previous elements or isotopes.  However, you may
+        use item assignment of the `~plasmapy.atomic.IonizationState`
+        instance to assign new ionic fractions one element or isotope
+        at a time.
+
+        Raises
+        ------
+        AtomicError
+            If the ionic fractions cannot be set.
+
+        TypeError
+            If ``inputs`` is not a `list`, `tuple`, or `dict` during
+            instantiation, or if ``inputs`` is not a `dict` when it is
+            being set.
+
+        """
+
+        # A potential problem is that using item assignment on the
+        # ionic_fractions attribute could cause the original attributes
+        # to be overwritten without checks being performed.  We might
+        # eventually want to create a new class or subclass of UserDict
+        # that goes through these checks.  In the meantime, we should
+        # make it clear to users to set ionic_fractions by using item
+        # assignment on the IonizationStates instance as a whole.  An
+        # example of the problem is `s = IonizationStates(["He"])` being
+        # followed by `s.ionic_fractions["He"] = 0.3`.
+
+        if hasattr(self, '_ionic_fractions'):
+            if not isinstance(inputs, dict):
+                raise TypeError(
+                    "Can only reset ionic_fractions with a dict if "
+                    "ionic_fractions has been set already.")
+            old_particles = set(self.base_particles)
+            new_particles = {particle_symbol(key) for key in inputs.keys()}
+            missing_particles = old_particles - new_particles
+            if missing_particles:
+                raise AtomicError(
+                    "Can only reset ionic fractions with a dict if "
+                    "the new base particles are a superset of the "
+                    "prior base particles.  To change ionic fractions "
+                    "for one base particle, use item assignment on the "
+                    "IonizationStates instance instead.")
+
+        if isinstance(inputs, dict):
+            original_keys = inputs.keys()
+            ionfrac_types = {type(inputs[key]) for key in original_keys}
+            inputs_have_quantities = u.Quantity in ionfrac_types
+
+            if inputs_have_quantities and len(ionfrac_types) != 1:
+                raise TypeError(
+                    "Ionic fraction information may only be inputted "
+                    "as a Quantity object if all ionic fractions are "
+                    "Quantity arrays with units of inverse volume.")
+
+            try:
+                particles = {key: Particle(key) for key in original_keys}
+            except (InvalidParticleError, TypeError) as exc:
+                raise AtomicError(
+                    "Unable to create IonizationStates instance "
+                    "because not all particles are valid.") from exc
+
+            # The particles whose ionization states are to be recorded
+            # should be elements or isotopes but not ions or neutrals.
+
+            for key in particles.keys():
+                is_element = particles[key].is_category('element')
+                has_charge_info = particles[key].is_category(any_of=["charged", "uncharged"])
+
+                if not is_element or has_charge_info:
+                    raise AtomicError(
+                        f"{key} is not an element or isotope without "
+                        f"charge information.")
+
+            # We are sorting the elements/isotopes by atomic number and
+            # mass number since we will often want to plot and analyze
+            # things and this is the most sensible order.
+
+            sorted_keys = sorted(original_keys, key=lambda k: (
+                particles[k].atomic_number,
+                particles[k].mass_number if particles[k].isotope else 0,
+            ))
+
+            _elements_and_isotopes = []
+            _particle_instances = []
+            new_ionic_fractions = {}
+
+            if inputs_have_quantities:
+                n_elems = {}
+
+            for key in sorted_keys:
+                new_key = particles[key].particle
+                _particle_instances.append(particles[key])
+                if new_key in _elements_and_isotopes:
+                    raise AtomicError("Repeated particles in IonizationStates.")
+
+                nstates_input = len(inputs[key])
+                nstates = particles[key].atomic_number + 1
+                if nstates != nstates_input:
+                    raise AtomicError(
+                        f"The ionic fractions array for {key} must "
+                        f"have a length of {nstates}.")
+
+                _elements_and_isotopes.append(new_key)
+                if inputs_have_quantities:
+                    try:
+                        number_densities = inputs[key].to(u.m ** -3)
+                        n_elem = np.sum(number_densities)
+                        new_ionic_fractions[new_key] = np.array(number_densities / n_elem)
+                        n_elems[key] = n_elem
+                    except u.UnitConversionError as exc:
+                        raise AtomicError("Units are not inverse volume.") from exc
+                elif isinstance(inputs[key], np.ndarray) and inputs[key].dtype.kind == 'f':
+                    new_ionic_fractions[particles[key].particle] = inputs[key]
+                else:
+                    try:
+                        new_ionic_fractions[particles[key].particle] = \
+                            np.array(inputs[key], dtype=np.float)
+                    except ValueError as exc:
+                        raise AtomicError(f"Inappropriate ionic fractions for {key}.") from exc
+
+            for key in _elements_and_isotopes:
+                fractions = new_ionic_fractions[key]
+                if not np.all(np.isnan(fractions)):
+                    if np.min(fractions) < 0 or np.max(fractions) > 1:
+                        raise AtomicError(f"Ionic fractions for {key} are not between 0 and 1.")
+                    if not np.isclose(np.sum(fractions), 1, atol=self.tol, rtol=0):
+                        raise AtomicError(f"Ionic fractions for {key} are not normalized to 1.")
+
+            # When the inputs provide the densities, the abundances must
+            # not have been provided because that would be redundant
+            # or contradictory information.  The number density scaling
+            # factor might or might not have been provided.  Have the
+            # number density scaling factor default to the total number
+            # of neutrals and ions across all elements and isotopes, if
+            # it was not provided.  Then go ahead and calculate the
+            # abundances based on that.  However, we need to be careful
+            # that the abundances are not overwritten during the
+            # instantiation of the class.
+
+            if inputs_have_quantities:
+                if np.isnan(self.n):
+                    new_n = 0 * u.m ** -3
+                    for key in _elements_and_isotopes:
+                        new_n += n_elems[key]
+                    self.n = new_n
+
+                new_abundances = {}
+                for key in _elements_and_isotopes:
+                    new_abundances[key] = np.float(n_elems[key] / self.n)
+
+                self._pars['abundances'] = new_abundances
+
+        elif isinstance(inputs, (list, tuple)):
+
+            try:
+                _particle_instances = [Particle(particle) for particle in inputs]
+            except (InvalidParticleError, TypeError) as exc:
+                raise AtomicError("Invalid inputs to IonizationStates.") from exc
+
+            _particle_instances.sort(
+                key=lambda p: (p.atomic_number, p.mass_number if p.isotope else 0)
+            )
+            _elements_and_isotopes = [particle.particle for particle in _particle_instances]
+            new_ionic_fractions = {
+                particle.particle: np.full(
+                    particle.atomic_number + 1,
+                    fill_value=np.nan,
+                    dtype=np.float64
+                ) for particle in _particle_instances
+            }
+        else:
+            raise TypeError("Incorrect inputs to set ionic_fractions.")
+
+        for i in range(1, len(_particle_instances)):
+            if _particle_instances[i - 1].element == _particle_instances[i].element:
+                if not _particle_instances[i - 1].isotope and _particle_instances[i].isotope:
+                    raise AtomicError("Cannot have an element and isotopes of that element.")
+
+        self._particle_instances = _particle_instances
+        self._base_particles = _elements_and_isotopes
+        self._ionic_fractions = new_ionic_fractions
+
+    def normalize(self) -> None:
+        """
+        Normalize the ionic fractions so that the sum for each element
+        equals one.
+        """
+        for particle in self.base_particles:
+            tot = np.sum(self.ionic_fractions[particle])
+            self.ionic_fractions[particle] = self.ionic_fractions[particle] / tot
+
+    def equilibrate(
+            self,
+            T_e: u.K = np.nan * u.K,
+            particles: str = 'all',
+            kappa: Real = np.inf):
+        """
+        Set the ionic fractions to collisional ionization equilibrium.
+        Not implemented.
+
+        The electron temperature used to calculate the new equilibrium
+        ionic fractions will be the argument ``T_e`` to this method if
+        given, and otherwise the attribute ``T_e`` if no electon
+        temperature is provided to this method.
+
+        Parameters
+        ----------
+        T_e: ~astropy.units.Quantity, optional
+            The electron temperature.
+
+        particles: `list`, `tuple`, or `str`, optional
+            The elements and isotopes to be equilibrated.  If
+            ``particles`` is ``'all'`` (default), then all
+            elements and isotopes will be equilibrated.
+
+        kappa: Real
+            The value of kappa for a kappa distribution for electrons.
+
+        """
+        raise NotImplementedError
+
+    @property
+    @u.quantity_input
+    def n_e(self) -> u.m ** -3:
+        """
+        Return the electron number density under the assumption of
+        quasineutrality.
+        """
+        number_densities = self.number_densities
+        n_e = 0.0 * u.m ** -3
+        for elem in self.base_particles:
+            atomic_numb = atomic_number(elem)
+            number_of_ionization_states = atomic_numb + 1
+            integer_charges = np.linspace(0, atomic_numb, number_of_ionization_states)
+            n_e += np.sum(number_densities[elem] * integer_charges)
+        return n_e
+
+    @property
+    @u.quantity_input
+    def n(self) -> u.m ** -3:
+        """Return the number density scaling factor."""
+        return self._pars['n']
+
+    @n.setter
+    @u.quantity_input
+    def n(self, n: u.m ** -3):
+        """Set the number density scaling factor."""
+        try:
+            n = n.to(u.m ** -3)
+        except u.UnitConversionError as exc:
+            raise AtomicError("Units cannot be converted to u.m ** -3.") from exc
+        except Exception as exc:
+            raise AtomicError(f"{n} is not a valid number density.") from exc
+        if n < 0 * u.m ** -3:
+            raise AtomicError("Number density cannot be negative.")
+        self._pars['n'] = n.to(u.m ** -3)
+
+    @property
+    def number_densities(self) -> Dict[str, u.Quantity]:
+        """
+        Return a `dict` containing the number densities for element or
+        isotope.
+        """
+        return {
+            elem: self.n * self.abundances[elem] * self.ionic_fractions[elem]
+            for elem in self.base_particles
+        }
+
+    @property
+    def abundances(self) -> Optional[Dict]:
+        """Return the elemental abundances."""
+        return self._pars['abundances']
+
+    @abundances.setter
+    def abundances(self, abundances_dict: Optional[Dict]):
+        """
+        Set the elemental (or isotopic) abundances.  The elements and
+        isotopes must be the same as or a superset of the elements whose
+        ionization states are being tracked.
+        """
+        if abundances_dict is None:
+            self._pars['abundances'] = {elem: np.nan for elem in self.base_particles}
+        elif not isinstance(abundances_dict, dict):
+            raise TypeError(
+                f"The abundances attribute must be a dict with "
+                f"elements or isotopes as keys and real numbers "
+                f"representing relative abundances as values.")
+        else:
+            old_keys = abundances_dict.keys()
+            try:
+                new_keys_dict = {particle_symbol(old_key): old_key for old_key in old_keys}
+            except Exception:
+                raise AtomicError(
+                    f"The key {repr(old_key)} in the abundances "
+                    f"dictionary is not a valid element or isotope.")
+
+            new_elements = new_keys_dict.keys()
+
+            old_elements_set = set(self.base_particles)
+            new_elements_set = set(new_elements)
+
+            if old_elements_set - new_elements_set:
+                raise AtomicError(
+                    f"The abundances of the following particles are "
+                    f"missing: {old_elements_set - new_elements_set}")
+
+            new_abundances_dict = {}
+
+            for element in new_elements:
+                inputted_abundance = abundances_dict[new_keys_dict[element]]
+                try:
+                    inputted_abundance = float(inputted_abundance)
+                except Exception:
+                    raise TypeError(
+                        f"The abundance for {element} was provided as"
+                        f"{inputted_abundance}, which cannot be "
+                        f"converted to a real number.") from None
+
+                if inputted_abundance < 0:
+                    raise AtomicError(f"The abundance of {element} is negative.")
+                new_abundances_dict[element] = inputted_abundance
+
+            self._pars['abundances'] = new_abundances_dict
+
+    @property
+    def log_abundances(self) -> Dict[str, Real]:
+        """
+        Return a `dict` with atomic or isotope symbols as keys and the
+        base 10 logarithms of the relative abundances as the
+        corresponding values.
+        """
+        log_abundances_dict = {}
+        for key in self.abundances.keys():
+            log_abundances_dict[key] = np.log10(self.abundances[key])
+        return log_abundances_dict
+
+    @log_abundances.setter
+    def log_abundances(self, value: Optional[Dict[str, Real]]):
+        """
+        Set the base 10 logarithm of the relative abundances.
+        """
+        if value is not None:
+            try:
+                new_abundances_input = {}
+                for key in value.keys():
+                    new_abundances_input[key] = 10 ** value[key]
+                self.abundances = new_abundances_input
+            except Exception:
+                raise AtomicError("Invalid log_abundances.") from None
+
+    @property
+    @u.quantity_input(equivalencies=u.temperature_energy())
+    def T_e(self) -> u.K:
+        """Return the electron temperature."""
+        return self._pars['T_e']
+
+    @T_e.setter
+    @u.quantity_input(equivalencies=u.temperature_energy())
+    def T_e(self, electron_temperature: u.K):
+        """Set the electron temperature."""
+        try:
+            temperature = electron_temperature.to(u.K, equivalencies=u.temperature_energy())
+        except (AttributeError, u.UnitsError):
+            raise AtomicError(
+                f"{electron_temperature} is not a valid temperature.") from None
+        if temperature < 0 * u.K:
+            raise AtomicError("The electron temperature cannot be negative.")
+        self._pars['T_e'] = temperature
+
+    @property
+    def kappa(self) -> np.real:
+        """
+        Return the kappa parameter for a kappa distribution function
+        for electrons.
+
+        The value of ``kappa`` must be greater than ``1.5`` in order to
+        have a valid distribution function.  If ``kappa`` equals
+        `~numpy.inf`, then the distribution function reduces to a
+        Maxwellian.
+
+        """
+        return self._pars['kappa']
+
+    @kappa.setter
+    def kappa(self, value: Real):
+        """
+        Set the kappa parameter for a kappa distribution function for
+        electrons.  The value must be between ``1.5`` and `~numpy.inf`.
+        """
+        kappa_errmsg = "kappa must be a real number greater than 1.5"
+        if not isinstance(value, Real):
+            raise TypeError(kappa_errmsg)
+        if value <= 1.5:
+            raise ValueError(kappa_errmsg)
+        self._pars['kappa'] = np.real(value)
+
+    @property
+    def base_particles(self) -> List[str]:
+        """
+        Return a list of the elements and isotopes whose ionization
+        states are being kept track of.
+        """
+        return self._base_particles
+
+    @property
+    def tol(self) -> np.real:
+        """Return the absolute tolerance for comparisons."""
+        return self._tol
+
+    @tol.setter
+    def tol(self, atol: Real):
+        """
+        Set the absolute tolerance for comparisons.
+        """
+        if not isinstance(atol, Real):
+            raise TypeError("The attribute tol must be a real number.")
+        if 0 <= atol <= 1.0:
+            self._tol = np.real(atol)
+        else:
+            raise ValueError("Need 0 <= tol <= 1.")
+
+    def info(self, minimum_ionic_fraction: Real = 0.01) -> None:
+        """
+        Print quicklook information for an
+        `~plasmapy.atomic.IonizationStates` instance.
+
+        Parameters
+        ----------
+        minimum_ionic_fraction: Real
+            If the ionic fraction for a particular ionization state is
+            below this level, then information for it will not be
+            printed.  Defaults to 0.01.
+
+        Examples
+        --------
+        >>> states = IonizationStates(
+        ...     {'H': [0.1, 0.9], 'He': [0.95, 0.05, 0.0]},
+        ...     T_e = 12000 * u.K,
+        ...     n = 3e9 * u.cm ** -3,
+        ...     abundances = {'H': 1.0, 'He': 0.1},
+        ...     kappa = 3.4,
+        ... )
+        >>> states.info()
+        IonizationStates instance for: H, He
+        ----------------------------------------------------------------
+        H  0+: 0.100    n_i = 3.00e+14 m**-3
+        H  1+: 0.900    n_i = 2.70e+15 m**-3
+        ----------------------------------------------------------------
+        He  0+: 0.950    n_i = 2.85e+14 m**-3
+        He  1+: 0.050    n_i = 1.50e+13 m**-3
+        ----------------------------------------------------------------
+        n_e = 2.71e+15 m**-3
+        T_e = 1.20e+04 K
+        kappa = 3.40
+        ----------------------------------------------------------------
+
+        """
+        separator_line = 64 * '-'
+
+        output = []
+
+        output.append(f"IonizationStates instance for: {', '.join(self.base_particles)}")
+
+        # Get the ionic symbol with the corresponding ionic fraction and
+        # number density (if available), but only for the most abundant
+        # ionization levels for each element.
+
+        for ionization_state in self:
+            states_info = ionization_state._get_states_info(minimum_ionic_fraction)
+            if len(states_info) > 0:
+                output += states_info
+                output[-1] += '\n' + separator_line
+
+        attributes = []
+        if np.isfinite(self.n_e):
+            attributes.append("n_e = " + "{:.2e}".format(self.n_e.value) + " m**-3")
+        if np.isfinite(self.T_e):
+            attributes.append("T_e = " + "{:.2e}".format(self.T_e.value) + " K")
+        if np.isfinite(self.kappa):
+            attributes.append("kappa = " + "{:.2f}".format(self.kappa))
+
+        if attributes:
+            attributes.append(separator_line)
+
+        output.append("\n".join(attributes))
+
+        if len(output) > 1:
+            output[0] += '\n' + separator_line
+            output_string = "\n".join(output)
+        else:
+            output_string = output[0]
+
+        print(output_string.strip('\n'))
diff --git a/plasmapy/atomic/particle_class.py b/plasmapy/atomic/particle_class.py
index bfa931d3..e89c4378 100644
--- a/plasmapy/atomic/particle_class.py
+++ b/plasmapy/atomic/particle_class.py
@@ -1,15 +1,17 @@
 """The Particle class."""
 
-import numpy as np
 import warnings
 from typing import (Union, Set, Tuple, List, Optional)
 import collections
-import plasmapy.utils.roman as roman
-import numbers
+from numbers import Integral, Real
 
+import numpy as np
 import astropy.units as u
 import astropy.constants as const
 
+import plasmapy.utils.roman as roman
+from plasmapy.atomic.elements import _Elements, _PeriodicTable
+from plasmapy.atomic.isotopes import _Isotopes
 from plasmapy.utils import (
     AtomicError,
     AtomicWarning,
@@ -21,17 +23,12 @@
     MissingAtomicDataError,
     MissingAtomicDataWarning,
 )
-
-from .parsing import (
+from plasmapy.atomic.parsing import (
     _dealias_particle_aliases,
     _parse_and_check_atomic_input,
     _invalid_particle_errmsg,
 )
-
-from .elements import _Elements, _PeriodicTable
-from .isotopes import _Isotopes
-
-from .special_particles import (
+from plasmapy.atomic.special_particles import (
     _Particles,
     ParticleZoo,
     _special_ion_masses,
@@ -39,7 +36,7 @@
 )
 
 __all__ = [
-    "Particle",
+    'Particle',
 ]
 
 _classification_categories = {
@@ -173,10 +170,8 @@ class Particle:
     >>> positron.particle
     'e+'
 
-    The `~plasmapy.atomic.Particle.atomic_symbol`,
-    `~plasmapy.atomic.Particle.isotope_symbol`, and
-    `~plasmapy.atomic.Particle.ionic_symbol` attributes return the
-    symbols for each of these different types of particles.
+    The ``element``, ``isotope``, and ``ionic_symbol`` attributes return
+    the symbols for each of these different types of particles.
 
     >>> proton.element
     'H'
@@ -280,19 +275,19 @@ class Particle:
 
     def __init__(
             self,
-            argument: Union[str, numbers.Integral],
-            mass_numb: numbers.Integral = None,
-            Z: numbers.Integral = None):
+            argument: Union[str, Integral],
+            mass_numb: Integral = None,
+            Z: Integral = None):
         """
         Instantiate a `~plasmapy.atomic.Particle` object and set private
         attributes.
         """
 
-        if not isinstance(argument, (numbers.Integral, np.integer, str, Particle)):
+        if not isinstance(argument, (Integral, np.integer, str, Particle)):
             raise TypeError(
-                "The first positional argument when creating a Particle "
-                "object must be either an integer, string, or another"
-                "Particle object.")
+                "The first positional argument when creating a "
+                "Particle object must be either an integer, string, or "
+                "another Particle object.")
 
         # If argument is a Particle instance, then we will construct a
         # new Particle instance for the same Particle (essentially a
@@ -301,10 +296,10 @@ def __init__(
         if isinstance(argument, Particle):
             argument = argument.particle
 
-        if mass_numb is not None and not isinstance(mass_numb, numbers.Integral):
+        if mass_numb is not None and not isinstance(mass_numb, Integral):
             raise TypeError("mass_numb is not an integer")
 
-        if Z is not None and not isinstance(Z, numbers.Integral):
+        if Z is not None and not isinstance(Z, Integral):
             raise TypeError("Z is not an integer.")
 
         self._attributes = collections.defaultdict(lambda: None)
@@ -527,7 +522,7 @@ def __ne__(self, other) -> bool:
         """
         Test whether or not two objects are different particles.
 
-        This method will return `False` if `other` is an identical
+        This method will return `False` if ``other`` is an identical
         `~plasmapy.atomic.Particle` instance or a `str` representing the
         same particle, and return `True` if ``other`` is a different
         `~plasmapy.atomic.Particle` or a `str` representing a different
@@ -542,12 +537,19 @@ def __ne__(self, other) -> bool:
         """
         return not self.__eq__(other)
 
+    def __hash__(self) -> int:
+        """
+        Allow use of `hash` so that a `~plasmapy.atomic.Particle`
+        instance may be used as a key in a `dict`.
+        """
+        return hash(self.__repr__())
+
     def __bool__(self):
         """
         Raise an `~plasmapy.utils.AtomicError` because Particle objects
         do not have a truth value.
         """
-        raise AtomicError("The truthiness of a Particle object is not defined.")
+        raise AtomicError("The truthiness of a Particle instance is not defined.")
 
     def __invert__(self):
         """
@@ -734,7 +736,7 @@ def isotope_name(self) -> str:
         return isotope_name
 
     @property
-    def integer_charge(self) -> numbers.Integral:
+    def integer_charge(self) -> Integral:
         """
         Return the particle's integer charge.
 
@@ -1032,7 +1034,7 @@ def binding_energy(self) -> u.Quantity:
         return nuclear_binding_energy.to(u.J)
 
     @property
-    def atomic_number(self) -> numbers.Integral:
+    def atomic_number(self) -> Integral:
         """
         Return the number of protons in an element, isotope, or ion.
 
@@ -1054,7 +1056,7 @@ def atomic_number(self) -> numbers.Integral:
         return self._attributes['atomic number']
 
     @property
-    def mass_number(self) -> numbers.Integral:
+    def mass_number(self) -> Integral:
         """
         Return the number of nucleons in an isotope.
 
@@ -1076,7 +1078,7 @@ def mass_number(self) -> numbers.Integral:
         return self._attributes['mass number']
 
     @property
-    def neutron_number(self) -> numbers.Integral:
+    def neutron_number(self) -> Integral:
         """
         Return the number of neutrons in an isotope or nucleon.
 
@@ -1103,7 +1105,7 @@ def neutron_number(self) -> numbers.Integral:
             raise InvalidIsotopeError(_category_errmsg(self, 'isotope'))
 
     @property
-    def electron_number(self) -> numbers.Integral:
+    def electron_number(self) -> Integral:
         """
         Return the number of electrons in an ion.
 
@@ -1161,7 +1163,7 @@ def isotopic_abundance(self) -> u.Quantity:
         return abundance
 
     @property
-    def baryon_number(self) -> numbers.Integral:
+    def baryon_number(self) -> Integral:
         """
         Return the number of baryons in a particle.
 
@@ -1185,7 +1187,7 @@ def baryon_number(self) -> numbers.Integral:
         return self._attributes['baryon number']
 
     @property
-    def lepton_number(self) -> numbers.Integral:
+    def lepton_number(self) -> Integral:
         """
         Return ``1`` for leptons, ``-1`` for antileptons, and ``0``
         otherwise.
@@ -1242,7 +1244,7 @@ def half_life(self) -> Union[u.Quantity, str]:
         return self._attributes['half-life']
 
     @property
-    def spin(self) -> Union[numbers.Integral, numbers.Real]:
+    def spin(self) -> Real:
         """
         Return the spin of the particle.
 
@@ -1290,7 +1292,8 @@ def periodic_table(self) -> collections.namedtuple:
 
     @property
     def categories(self) -> Set[str]:
-        """Return the particle's categories.
+        """
+        Return the particle's categories.
 
         Examples
         --------
@@ -1317,9 +1320,9 @@ def is_category(self,
 
         Required categories may be entered as positional arguments,
         including as a `list`, `set`, or `tuple` of required categories.
-        These may also be included using the `require` keyword argument.
-        This method will return `False` if the particle is not in all of
-        the required categories.
+        These may also be included using the ``require`` keyword
+        argument.  This method will return `False` if the particle is
+        not in all of the required categories.
 
         If categories are inputted using the ``any_of`` keyword
         argument, then this method will return `False` if the particle
@@ -1410,8 +1413,7 @@ def is_electron(self) -> bool:
     @property
     def is_ion(self) -> bool:
         """
-        Return `True` if the particle is an ion, and `False`
-        otherwise.
+        Return `True` if the particle is an ion, and `False` otherwise.
 
         Examples
         --------
@@ -1425,7 +1427,7 @@ def is_ion(self) -> bool:
         """
         return self.is_category('ion')
 
-    def ionize(self, n: numbers.Integral = 1, inplace: bool = False):
+    def ionize(self, n: Integral = 1, inplace: bool = False):
         """
         Create a new `~plasmapy.atomic.Particle` instance corresponding
         to the current `~plasmapy.atomic.Particle` after being ionized
@@ -1494,7 +1496,7 @@ def ionize(self, n: numbers.Integral = 1, inplace: bool = False):
             raise InvalidIonError(
                 f"The particle {self.particle} is already fully "
                 f"ionized and cannot be ionized further.")
-        if not isinstance(n, numbers.Integral):
+        if not isinstance(n, Integral):
             raise TypeError("n must be a positive integer.")
         if n <= 0:
             raise ValueError("n must be a positive number.")
@@ -1507,7 +1509,7 @@ def ionize(self, n: numbers.Integral = 1, inplace: bool = False):
         else:
             return Particle(base_particle, Z=new_integer_charge)
 
-    def recombine(self, n: numbers.Integral = 1, inplace=False):
+    def recombine(self, n: Integral = 1, inplace=False):
         """
         Create a new `~plasmapy.atomic.Particle` instance corresponding
         to the current `~plasmapy.atomic.Particle` after undergoing
@@ -1537,7 +1539,7 @@ def recombine(self, n: numbers.Integral = 1, inplace=False):
         particle : ~plasmapy.atomic.Particle
             A new `~plasmapy.atomic.Particle` object that has undergone
             recombination ``n`` times relative to the original
-            `~plasmapy.atomic.Particle`.  If `inplace` is `False`,
+            `~plasmapy.atomic.Particle`.  If ``inplace`` is `False`,
             instead return `None`.
 
         Raises
@@ -1571,7 +1573,7 @@ def recombine(self, n: numbers.Integral = 1, inplace=False):
             raise ChargeError(
                 f"{self.particle} cannot undergo recombination because "
                 f"its charge is not specified.")
-        if not isinstance(n, numbers.Integral):
+        if not isinstance(n, Integral):
             raise TypeError("n must be a positive integer.")
         if n <= 0:
             raise ValueError("n must be a positive number.")
diff --git a/plasmapy/atomic/particle_input.py b/plasmapy/atomic/particle_input.py
index e50d4172..f7f5bf10 100644
--- a/plasmapy/atomic/particle_input.py
+++ b/plasmapy/atomic/particle_input.py
@@ -11,12 +11,14 @@
 
 from .particle_class import Particle
 
-from plasmapy.utils import (AtomicError,
-                     InvalidParticleError,
-                     InvalidElementError,
-                     InvalidIonError,
-                     InvalidIsotopeError,
-                     ChargeError)
+from plasmapy.utils import (
+    AtomicError,
+    InvalidParticleError,
+    InvalidElementError,
+    InvalidIonError,
+    InvalidIsotopeError,
+    ChargeError,
+)
 
 __all__ = [
     "particle_input",
diff --git a/plasmapy/classes/sources/openpmd_hdf5.py b/plasmapy/classes/sources/openpmd_hdf5.py
index 8b4c602c..bf53b121 100644
--- a/plasmapy/classes/sources/openpmd_hdf5.py
+++ b/plasmapy/classes/sources/openpmd_hdf5.py
@@ -1,6 +1,10 @@
-import h5py
 import numpy as np
 import astropy.units as u
+try:
+    import h5py
+except (ImportError, ModuleNotFoundError) as e:
+    from plasmapy.optional_deps import h5py_import_error
+    raise ImportError(h5py_import_error) from e
 
 from plasmapy.classes import GenericPlasma
 from plasmapy.utils import DataStandardError
diff --git a/plasmapy/diagnostics/langmuir.py b/plasmapy/diagnostics/langmuir.py
index f645d37d..5d8fff7b 100644
--- a/plasmapy/diagnostics/langmuir.py
+++ b/plasmapy/diagnostics/langmuir.py
@@ -6,7 +6,11 @@
 import plasmapy.constants as const
 import copy
 from astropy.visualization import quantity_support
-import matplotlib.pyplot as plt
+try:
+    import matplotlib.pyplot as plt
+except (ImportError, ModuleNotFoundError) as e:
+    from plasmapy.optional_deps import mpl_import_error
+    raise mpl_import_error from e
 from scipy.optimize import curve_fit
 from plasmapy.atomic import Particle
 
diff --git a/plasmapy/mathematics/mathematics.py b/plasmapy/mathematics/mathematics.py
index 7474d6a7..6aa17a52 100644
--- a/plasmapy/mathematics/mathematics.py
+++ b/plasmapy/mathematics/mathematics.py
@@ -3,7 +3,11 @@
 import numbers
 import numpy as np
 from astropy import units as u
-from mpmath import polylog
+try:
+    from mpmath import polylog
+except (ImportError, ModuleNotFoundError) as e:
+    from plasmapy.optional_deps import mpmath_import_error
+    raise mpmath_import_error from e
 from scipy.special import wofz as Faddeeva_function
 from typing import Union
 
diff --git a/plasmapy/optional_deps.py b/plasmapy/optional_deps.py
new file mode 100644
index 00000000..1a360b32
--- /dev/null
+++ b/plasmapy/optional_deps.py
@@ -0,0 +1,61 @@
+""" Useful error messages for optional dependencies that aren't found. """
+from typing import Optional
+
+
+def _optional_import_error_template(pkgname: str,
+                                    url: str,
+                                    library: Optional[str] = None,
+                                    conda_channel: Optional[str] = None,
+                                    ) -> ImportError:
+    """
+    Simple template to prevent text duplication.
+
+    Parameters
+    ----------
+    pkgname : str
+        Package name on pip or conda
+    url : str
+        Link to install page for given package
+    library : str (optional)
+        Full package name
+    conda_channel: str (optional)
+        set to "conda-forge" to add -c conda-forge to the conda install line
+
+    Returns
+    -------
+    ImportError
+
+    """
+    library = pkgname if library is None else library
+    conda_channel = f"-c {conda_channel} " if conda_channel is not None else ""
+
+    template = f"""
+    {library} could not be found. Try either
+
+        conda install {conda_channel}{pkgname}
+        
+    on Anaconda environments or
+        
+        pip install {pkgname}
+
+    in general. In case of trouble refer to
+    
+        {url} 
+    
+    (link active as of 2018.10.31 - please report dead links on GitHub!)"""
+
+    return ImportError(template)
+
+
+h5py_import_error = _optional_import_error_template("h5py",
+                              "http://docs.h5py.org/en/latest/build.html")
+
+mpl_import_error = _optional_import_error_template("matplotlib",
+                             "https://matplotlib.org/users/installing.html")
+
+mpmath_import_error = _optional_import_error_template("mpmath",
+                                "http://mpmath.org/doc/current/setup.html#download-and-installation")
+
+lmfit_import_error = _optional_import_error_template("lmfit",
+                               "https://lmfit.github.io/lmfit-py/installation.html")
+
diff --git a/plasmapy/physics/parameters.py b/plasmapy/physics/parameters.py
index 86f78b29..07eb2d33 100644
--- a/plasmapy/physics/parameters.py
+++ b/plasmapy/physics/parameters.py
@@ -1154,28 +1154,28 @@ def Debye_number(T_e: u.K, n_e: u.m**-3):
     return N_D.to(u.dimensionless_unscaled)
 
 
+
 @utils.check_quantity(
     n={'units': u.m ** -3, 'can_be_negative': False}
     )
-def inertial_length(n: u.m**-3, particle='e-'):
-    r"""Calculate the particle inertial length. At this length, the Hall effect
-    becomes important.
+@atomic.particle_input(require='charged')
+def inertial_length(n: u.m**-3, particle: atomic.Particle):
+    r"""
+    Calculate a charged particle's inertial length.
 
     Parameters
     ----------
     n : ~astropy.units.Quantity
-        Particle number density in units convertible to m**-3.
+        Particle number density in units convertible to m ** -3.
 
     particle : str, optional
-        Representation of the particle species (e.g., 'p' for protons, 'D+'
-        for deuterium, or 'He-4 +1' for singly ionized helium-4),
-        which defaults to electrons.  If no charge state information is
-        provided, then the particles are assumed to be singly charged.
+        Representation of the particle species (e.g., 'p+' for protons,
+        'D+' for deuterium, or 'He-4 +1' for singly ionized helium-4).
 
     Returns
     -------
     d : ~astropy.units.Quantity
-        Particles inertial length in meters.
+        The particle's inertial length in meters.
 
     Raises
     ------
@@ -1191,22 +1191,27 @@ def inertial_length(n: u.m**-3, particle='e-'):
     Warns
     -----
     ~astropy.units.UnitsWarning
-        If units are not provided, SI units are assumed
+        If units are not provided and SI units are assumed.
 
     Notes
     -----
-    The particle inertial length is also known as an particle skin depth and is
-    given by:
+    The inertial length of a particle of species :math:`s` is given by
 
     .. math::
-        d = \frac{c}{\omega_{pi}}
+        d = \frac{c}{\omega_{ps}}
+
+    The inertial length is the characteristic length scale for a
+    particle to be accelerated in a plasma.  The Hall effect becomes
+    important on length scales shorter than the ion inertial length.
+
+    The inertial length is also known as the skin depth.
 
     Example
     -------
     >>> from astropy import units as u
-    >>> inertial_length(5*u.m**-3, particle='He+')
+    >>> inertial_length(5 * u.m ** -3, 'He+')
     <Quantity 2.02985802e+08 m>
-    >>> inertial_length(5*u.m**-3)
+    >>> inertial_length(5 * u.m ** -3, 'e-')
     <Quantity 2376534.75601976 m>
 
     """
diff --git a/plasmapy/physics/quantum.py b/plasmapy/physics/quantum.py
index 0cc7ef2e..a9a85f7b 100644
--- a/plasmapy/physics/quantum.py
+++ b/plasmapy/physics/quantum.py
@@ -6,6 +6,11 @@
 # python modules
 import numpy as np
 from astropy import units as u
+try:
+    from lmfit import minimize, Parameters
+except (ImportError, ModuleNotFoundError) as e:
+    from plasmapy.optional_deps import lmfit_import_error
+    raise lmfit_import_error from e
 
 # plasmapy modules
 from plasmapy import atomic, utils, mathematics
@@ -447,7 +452,6 @@ def chemical_potential(n_e: u.m ** -3, T: u.K):
     <Quantity 2.00039985e-12>
 
     """
-    from lmfit import minimize, Parameters
     # deBroglie wavelength
     lambdaDB = thermal_deBroglie_wavelength(T)
     # degeneracy parameter
diff --git a/setup.cfg b/setup.cfg
index f7316236..ef355f0c 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -9,7 +9,15 @@ edit_on_github = True
 github_project = PlasmaPy/PlasmaPy
 # install_requires should be formatted as a comma-separated list, e.g.:
 # install_requires = astropy, scipy, matplotlib
-install_requires = numpy, scipy, astropy, cython, lmfit, matplotlib, colorama, mpmath, h5py
+install_requires = numpy, scipy, astropy, cython, colorama
+
+extra_requires = classes, diagnostics, tests, math, physics
+classes_requires = h5py, matplotlib
+diagnostics_requires = matplotlib
+tests_requires = pytest, pytest_astropy
+math_requires = mpmath
+physics_requires = lmfit
+
 # version should be PEP386 compatible (http://www.python.org/dev/peps/pep-0386)
 version = 0.1.1.dev
 
diff --git a/setup.py b/setup.py
index 46c9ca79..25b6eb9d 100644
--- a/setup.py
+++ b/setup.py
@@ -4,6 +4,7 @@
 import glob
 import os
 import sys
+import itertools
 
 # Enforce Python version check - this is the same check as in __init__.py but
 # this one has to happen before importing ah_bootstrap.
@@ -127,6 +128,18 @@
 
 setup_requires = ['numpy']
 
+extra_tags = [m.strip() for m in metadata.get("extra_requires", "").split(',')] 
+if extra_tags:
+     # TODO: the line below will prove useful once #576 is in play
+     extras_require = {tag: [m.strip() for m in metadata[f"{tag}_requires"].split(',')] 
+                       for tag in extra_tags}
+
+    # but for now we'll limit ourselves to pip install plasmapy[optional]
+     optionals_require = list(itertools.chain.from_iterable(extras_require.values()))
+     extras_require = dict(optional=optionals_require)
+else:
+     extras_require = None 
+
 # Make sure to have the packages needed for building PlasmaPy, but do not require them
 # when installing from an sdist as the c files are included there.
 if not os.path.exists(os.path.join(os.path.dirname(__file__), 'PKG-INFO')):
@@ -172,6 +185,7 @@
       include_package_data=True,
       entry_points=entry_points,
       python_requires='>={}'.format("3.6"),
-      tests_require=["pytest", "pytest-astropy"],
+      tests_require=extras_require.get("optional", ""),
+      extras_require=extras_require,
       **package_info
 )


Test patch:
diff --git a/plasmapy/atomic/tests/test_ionization_state.py b/plasmapy/atomic/tests/test_ionization_state.py
new file mode 100644
index 00000000..7f935591
--- /dev/null
+++ b/plasmapy/atomic/tests/test_ionization_state.py
@@ -0,0 +1,553 @@
+import pytest
+import collections
+
+import numpy as np
+import astropy.units as u
+
+from plasmapy.atomic.ionization_state import IonizationState
+from plasmapy.utils import AtomicError, InvalidIsotopeError, run_test
+from plasmapy.atomic import (
+    atomic_number,
+    atomic_symbol,
+    particle_symbol,
+    isotope_symbol,
+    Particle,
+)
+
+test_cases = {
+    'Li': {
+        'particle': 'Li',
+        'ionic_fractions': np.array([0.4, 0.3, 0.2, 0.1]),
+        'tol': 1e-15,
+    },
+
+    'Li ground state': {
+        'particle': 'Li',
+        'ionic_fractions': np.array([1, 0, 0, 0], dtype=np.int64),
+        'tol': 1e-15,
+    },
+
+    'H': {
+        'particle': 'H',
+        'ionic_fractions': [0.6, 0.4],
+        'tol': 1e-8,
+    },
+
+    'H acceptable error': {
+        'particle': 'H',
+        'ionic_fractions': [0.6, 0.400_000_001],
+        'tol': 1e-8,
+    },
+
+    'D': {
+        'particle': 'deuterium',
+        'ionic_fractions': [0.7, 0.3],
+        'tol': 1e-15,
+        'n_elem': 3e14 * u.m ** -3,
+    },
+
+    'He': {
+        'particle': 'He',
+        'ionic_fractions': [0.5, 0.3, 0.2],
+        'n_elem': 1e20 * u.m ** -3,
+    },
+
+    'number densities': {
+        'particle': 'T',
+        'ionic_fractions': np.array([1e4, 1]) * u.cm ** -3,
+    }
+
+}
+
+test_names = test_cases.keys()
+
+
+class Test_IonizationState:
+    """Test instances of IonizationState."""
+
+    @classmethod
+    def setup_class(cls):
+        cls.instances = {}
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_instantiation(self, test_name):
+        """
+        Test that each IonizationState test case can be instantiated.
+        """
+        try:
+            self.instances[test_name] = IonizationState(**test_cases[test_name])
+        except Exception:
+            pytest.fail(f"Unable to create IonizationState instance for test case {test_name}.")
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_integer_charges(self, test_name):
+        """
+        Test that an `IonizationState` instance has the correct integer
+        charges.
+        """
+        instance = self.instances[test_name]
+        expected_integer_charges = np.arange(instance.atomic_number + 1)
+        assert np.allclose(instance.integer_charges, expected_integer_charges)
+
+    @pytest.mark.parametrize(
+        'test_name',
+        [name for name in test_names if 'ionic_fractions' in test_cases[name].keys()],
+    )
+    def test_ionic_fractions(self, test_name):
+        """
+        Test that each `IonizationState` instance has the expected
+        ionic fractions.
+        """
+        instance = self.instances[test_name]
+        inputted_fractions = test_cases[test_name]['ionic_fractions']
+        if isinstance(inputted_fractions, u.Quantity):
+            inputted_fractions = inputted_fractions.to(u.m ** -3)
+            inputted_fractions = (inputted_fractions / inputted_fractions.sum()).value
+        if not np.allclose(instance.ionic_fractions, inputted_fractions):
+            pytest.fail(f"Mismatch in ionic fractions for test {test_name}.")
+
+    def test_equal_to_itself(self):
+        """
+        Test that `IonizationState.__eq__` returns `True for two identical
+        `IonizationState` instances.
+        """
+        assert self.instances['Li'] == self.instances['Li'], \
+            "Identical IonizationState instances are not equal."
+
+    def test_equal_to_within_tolerance(self):
+        """
+        Test that `IonizationState.__eq__` returns `True` for two
+        `IonizationState` instances that differ within the inputted
+        tolerance.
+        """
+        assert self.instances['H'] == self.instances['H acceptable error'], \
+            ("Two IonizationState instances that are approximately the "
+             "same to within the tolerance are not testing as equal.")
+
+    def test_inequality(self):
+        """
+        Test that instances with different ionic fractions are not equal
+        to each other.
+        """
+        assert self.instances['Li ground state'] != self.instances['Li'], \
+            "Different IonizationState instances are equal."
+
+    def test_equality_exception(self):
+        """
+        Test that comparisons of `IonizationState` instances for
+        different elements fail.
+        """
+        with pytest.raises(AtomicError):
+            self.instances['Li'] == self.instances['H']
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_iteration(self, test_name: str):
+        """Test that `IonizationState` instances iterate impeccably."""
+        try:
+            states = [state for state in self.instances[test_name]]
+        except Exception:
+            pytest.fail(f"Unable to perform iteration for {test_name}.")
+
+        try:
+            integer_charges = [state.integer_charge for state in states]
+            ionic_fractions = np.array([state.ionic_fraction for state in states])
+            ionic_symbols = [state.ionic_symbol for state in states]
+        except Exception:
+            pytest.fail(f"An attribute may be misnamed or missing ({test_name}).")
+
+        try:
+            base_symbol = isotope_symbol(ionic_symbols[0])
+        except InvalidIsotopeError:
+            base_symbol = atomic_symbol(ionic_symbols[0])
+        finally:
+            atomic_numb = atomic_number(ionic_symbols[1])
+
+        errors = []
+
+        expected_charges = np.arange(atomic_numb + 1)
+        if not np.all(integer_charges == expected_charges):
+            errors.append(
+                f"The resulting integer charges are {integer_charges}, "
+                f"which are not equal to the expected integer charges, "
+                f"which are {expected_charges}.")
+
+        expected_fracs = test_cases[test_name]['ionic_fractions']
+        if isinstance(expected_fracs, u.Quantity):
+            expected_fracs = (expected_fracs / expected_fracs.sum()).value
+
+        if not np.allclose(ionic_fractions, expected_fracs):
+            errors.append(
+                f"The resulting ionic fractions are {ionic_fractions}, "
+                f"which are not equal to the expected ionic fractions "
+                f"of {expected_fracs}.")
+
+        expected_particles = [Particle(base_symbol, Z=charge) for charge in integer_charges]
+        expected_symbols = [particle.ionic_symbol for particle in expected_particles]
+        if not ionic_symbols == expected_symbols:
+            errors.append(
+                f"The resulting ionic symbols are {ionic_symbols}, "
+                f"which are not equal to the expected ionic symbols of "
+                f"{expected_symbols}.")
+
+        if errors:
+            errors.insert(0, (
+                f"The test of IonizationState named '{test_name}' has "
+                f"resulted in the following errors when attempting to "
+                f"iterate."))
+            errmsg = " ".join(errors)
+            pytest.fail(errmsg)
+
+    def test_slicing_error(self):
+        """
+        Test that an IonizationState instance cannot be sliced.
+        """
+        with pytest.raises(TypeError):
+            self.instances['Li'][1:3]
+
+    @pytest.mark.parametrize('index', [-1, 4, 'Li'])
+    def test_indexing_error(self, index):
+        """
+        Test that an `IonizationState` instance cannot be indexed
+        outside of the bounds of allowed integer charges.
+        """
+        with pytest.raises(AtomicError):
+            self.instances['Li'][index]
+
+    def test_normalization(self):
+        """
+        Test that `_is_normalized` returns `False` when there is an
+        error greater than the tolerance, and `True` after normalizing.
+        """
+        H = self.instances['H acceptable error']
+        assert not H._is_normalized(tol=1e-15)
+        H.normalize()
+        assert H._is_normalized(tol=1e-15)
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_identifications(self, test_name):
+        """
+        Test that the identification attributes for test
+        `IonizationState` instances match the expected values from the
+        `Particle` instance.
+        """
+
+        Identifications = collections.namedtuple(
+            "Identifications",
+            ["element", "isotope", "atomic_number"],
+        )
+
+        expected_identifications = Identifications(
+            self.instances[test_name].element,
+            self.instances[test_name].isotope,
+            self.instances[test_name].atomic_number,
+        )
+
+        expected_element = self.instances[test_name]._particle_instance.element
+        expected_isotope = self.instances[test_name]._particle_instance.isotope
+        expected_atomic_number = self.instances[test_name]._particle_instance.atomic_number
+
+        resulting_identifications = Identifications(
+            expected_element,
+            expected_isotope,
+            expected_atomic_number,
+        )
+
+        assert resulting_identifications == expected_identifications, (
+            f"For IonizationState test {test_name}, the resulting "
+            f"identifications of {resulting_identifications} differ "
+            f"from the expected identifications of "
+            f"{expected_identifications}."
+        )
+
+    @pytest.mark.parametrize('tol', [-1e-16, 1.0000001])
+    def test_invalid_tolerances(self, tol):
+        """Test that invalid tolerances raise appropriate errors."""
+        test_name = "Li"
+        instance = self.instances[test_name]
+        with pytest.raises(ValueError):
+            instance.tol = tol
+
+    @pytest.mark.parametrize('test_name', test_cases.keys())
+    def test_particle_instances(self, test_name):
+        """
+        Test that `IonizationState` returns the correct `Particle`
+        instances.
+        """
+        instance = self.instances[test_name]
+        atom = instance.base_particle
+        nstates = instance.atomic_number + 1
+        expected_particles = [Particle(atom, Z=Z) for Z in range(nstates)]
+        assert expected_particles == instance._particle_instances, (
+            f"The expected Particle instances of {expected_particles} "
+            f"are not all equal to the IonizationState particles of "
+            f"{instance._particle_instances} for test {test_name}.")
+
+    @pytest.mark.parametrize(
+        'test_name',
+        [name for name in test_names if 'n_elem' in test_cases[name].keys()])
+    def test_electron_density_from_n_elem_ionic_fractions(self, test_name):
+        instance = self.instances[test_name]
+        n_elem = test_cases[test_name]['n_elem']
+        ionic_fractions = test_cases[test_name]['ionic_fractions']
+        assert instance.n_elem == n_elem, \
+            f"n_elem is not being stored correctly for test {test_name}"
+        assert np.isclose(
+            instance.n_e,
+            np.sum(n_elem * ionic_fractions * np.arange(instance.atomic_number + 1)),
+            rtol=1e-12, atol=0 * u.m ** -3), \
+            "n_e is not the expected value."
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_getitem(self, test_name):
+        """
+        Test that `IonizationState.__getitem__` returns the same value
+        when using equivalent keys (integer charge, particle symbol, and
+        `Particle` instance).
+
+        For example, if we create
+
+        >>> He_states = IonizationState('He', [0.2, 0.3, 0.5])
+
+        then this checks to make sure that `He_states[2]`,
+        `He_states['He 2+']`, and `He_states[Particle('He 2+')]` all
+        return the same result.
+
+        """
+        instance = self.instances[test_name]
+        particle_name = instance.base_particle
+
+        integer_charges = np.arange(instance.atomic_number + 1)
+        symbols = [particle_symbol(particle_name, Z=Z) for Z in integer_charges]
+        particles = instance._particle_instances
+
+        errors = []
+
+        # In the following loop, instance[key] will return a namedtuple
+        # or class which may contain Quantity objects with values of
+        # numpy.nan.  Because of the difficulty of comparing nans in
+        # these objects, we compare the string representations instead
+        # (see Astropy issue #7901 on GitHub).
+
+        for keys in zip(integer_charges, symbols, particles):
+            set_of_str_values = {str(instance[key]) for key in keys}
+            if len(set_of_str_values) != 1:
+                errors.append(
+                    f'\n\n'
+                    f'The following keys in test {test_name} did not '
+                    f'produce identical outputs as required: {keys}. '
+                    f'The set containing string representations of'
+                    f'the values is:\n\n{set_of_str_values}')
+
+        if errors:
+            pytest.fail(str.join('', errors))
+
+    @pytest.mark.parametrize('attr', ['integer_charge', 'ionic_fraction', 'ionic_symbol'])
+    def test_State_attrs(self, attr):
+        """
+        Test that an IonizationState instance returns something with the
+        correct attributes (be it a `collections.namedtuple` or a
+        `class`).
+        """
+        test_name = 'He'
+        state = self.instances[test_name][1]
+        assert hasattr(state, attr)
+
+    def test_State_equality_and_getitem(self):
+        test_name = 'He'
+        instance = self.instances[test_name]
+        charge = 2
+        symbol = 'He 2+'
+        result_from_charge = instance[charge]
+        result_from_symbol = instance[symbol]
+        assert result_from_charge == result_from_symbol
+
+
+IE = collections.namedtuple("IE", ["inputs", "expected_exception"])
+
+tests_for_exceptions = {
+    'too few nstates': IE({'particle': 'H', 'ionic_fractions': [1.0]}, AtomicError),
+    'too many nstates': IE({'particle': 'H', 'ionic_fractions': [1, 0, 0, 0]}, AtomicError),
+    'ionic fraction < 0': IE({'particle': 'He', 'ionic_fractions': [-0.1, 0.1, 1]}, AtomicError),
+    'ionic fraction > 1': IE({'particle': 'He', 'ionic_fractions': [1.1, 0.0, 0.0]}, AtomicError),
+
+    'invalid ionic fraction': IE({
+        'particle': 'He', 'ionic_fractions': [1.0, 0.0, 'a']
+    }, AtomicError),
+
+    'bad n_elem units': IE({
+        'particle': 'H', 'ionic_fractions': [0, 1], 'n_elem': 3 * u.m ** 3
+    }, u.UnitConversionError),
+
+    'bad T_e units': IE({
+        'particle': 'H', 'ionic_fractions': [0, 1], 'T_e': 1 * u.m
+    }, u.UnitConversionError),
+
+    'negative n_elem': IE({
+        'particle': 'He', 'ionic_fractions': [1.0, 0.0, 0.0], 'n_elem': -1 * u.m ** -3
+    }, AtomicError),
+
+    'negative T_e': IE({
+        'particle': 'He', 'ionic_fractions': [1.0, 0.0, 0.0], 'T_e': -1 * u.K
+    }, AtomicError),
+
+    'redundant ndens': IE({
+        'particle': 'H', 'ionic_fractions': np.array([3, 4]) * u.m ** -3, 'n_elem': 4 * u.m ** -3,
+    }, AtomicError),
+
+}
+
+
+@pytest.mark.parametrize('test', tests_for_exceptions.keys())
+def test_IonizationState_exceptions(test):
+    """
+    Test that appropriate exceptions are raised for inappropriate inputs
+    to `IonizationState`.
+    """
+    run_test(
+        IonizationState,
+        kwargs=tests_for_exceptions[test].inputs,
+        expected_outcome=tests_for_exceptions[test].expected_exception,
+    )
+
+
+kwargs = {
+    'particle': 'He-4',
+    'ionic_fractions': [0.2, 0.3, 0.5],
+    'T_e': 5.0 * u.kK,
+    'tol': 2e-14,
+    'n_elem': 1e13 * u.cm ** -3
+}
+
+expected_properties = {
+    'T_e': 5000.0 * u.K,
+    'tol': 2e-14,
+    'isotope': 'He-4',
+    'element': 'He',
+    'atomic_number': 2,
+    'Z_mean': 1.3,
+    'Z_rms': 1.51657508881031,
+    'n_e': 1.3e19 * u.m ** -3,
+    'n_elem': 1e19 * u.m ** -3,
+    'integer_charges': [0, 1, 2],
+    'ionic_fractions': np.array([0.2, 0.3, 0.5]),
+    'ionic_symbols': ['He-4 0+', 'He-4 1+', 'He-4 2+'],
+    '_is_normalized()': True,
+    'number_densities': np.array([2e18, 3e18, 5e18]) * u.m ** -3,
+    'tol': 2e-14,
+    '__str__()': "<IonizationState instance for He-4>",
+}
+
+instance = IonizationState(**kwargs)
+
+
+@pytest.mark.parametrize('key', expected_properties.keys())
+def test_IonizationState_attributes(key):
+    """
+    Test a specific case that the `IonizationState` attributes are
+    working as expected.
+    """
+    expected = expected_properties[key]
+    actual = eval(f'instance.{key}')
+
+    if isinstance(expected, u.Quantity):
+        assert expected.unit == actual.unit, f"Unit mismatch for IonizationState.{key}"
+        assert np.allclose(expected, actual, atol=1e-15 * expected.unit), \
+            f"Quantity.value mismatch for IonizationState.{key}"
+    else:
+        try:
+            assert expected == actual
+        except ValueError:
+            assert np.allclose(expected, actual)
+
+
+def test_nans():
+    """
+    Test that when no ionic fractions or temperature are inputted,
+    the result is an array full of `~numpy.nan` of the right size.
+    """
+    element = 'He'
+    nstates = atomic_number(element) + 1
+    instance = IonizationState(element)
+    assert len(instance.ionic_fractions) == nstates, \
+        f"Incorrect number of ionization states for {element}"
+    assert np.all([np.isnan(instance.ionic_fractions)]), (
+        f"The ionic fractions for IonizationState are not defaulting " 
+        f"to numpy.nan when not set by user.")
+
+
+def test_setting_ionic_fractions():
+    instance = IonizationState('He')
+    new_ionic_fractions = [0.2, 0.5, 0.3]
+    instance.ionic_fractions = new_ionic_fractions
+    assert np.allclose(instance.ionic_fractions, new_ionic_fractions)
+
+
+class Test_IonizationStateNumberDensitiesSetter:
+    """Test that setting IonizationState.number_densities works correctly."""
+
+    def setup_class(self):
+        self.element = 'H'
+        self.valid_number_densities = u.Quantity([0.1, 0.2], unit=u.m**-3)
+        self.expected_n_elem = np.sum(self.valid_number_densities)
+        self.expected_ionic_fractions = self.valid_number_densities / self.expected_n_elem
+        try:
+            self.instance = IonizationState(self.element)
+        except Exception:
+            pytest.fail("Unable to instantiate IonizationState with no ionic fractions.")
+
+    def test_setting_number_densities(self):
+        try:
+            self.instance.number_densities = self.valid_number_densities
+        except Exception:
+            pytest.fail(
+                f"Unable to set number densities of {self.element} to "
+                f"{self.valid_number_densities}.")
+
+        assert u.quantity.allclose(
+            self.instance.number_densities,
+            self.valid_number_densities), (
+            f"The number densities of {self.element} were set to "
+            f"{self.instance.number_densities} instead of the expceted "
+            f"value of {self.valid_number_densities}.")
+
+    def test_ionic_fractions(self):
+        assert np.allclose(
+            self.instance.ionic_fractions,
+            self.expected_ionic_fractions), (
+            "The IonizationState.ionic_fractions attribute was not set "
+            "correctly after the number densities were set.")
+
+    def test_n_elem(self):
+        assert u.quantity.allclose(
+            self.instance.n_elem,
+            self.expected_n_elem), (
+            "IonizationState.n_elem not set correctly after "
+            "number_densities was set.")
+
+    def test_n_e(self):
+        assert u.quantity.allclose(
+            self.instance.n_e,
+            self.valid_number_densities[1]), (
+            "IonizationState.n_e not set correctly after "
+            "number_densities was set.")
+
+    def test_that_negative_density_raises_error(self):
+        with pytest.raises(AtomicError, message="cannot be negative"):
+            self.instance.number_densities = u.Quantity([-0.1, 0.2], unit=u.m**-3)
+
+    def test_incorrect_number_of_charge_states_error(self):
+        with pytest.raises(AtomicError, message="Incorrect number of charge states"):
+            self.instance.number_densities = u.Quantity([0.1, 0.2, 0.3], unit=u.m**-3)
+
+    def test_incorrect_units_error(self):
+        with pytest.raises(u.UnitsError):
+            self.instance.number_densities = u.Quantity([0.1, 0.2], unit=u.kg)
+
+    # The following two tests are not related to setting the
+    # number_densities attribute, but are helpful to test anyway.
+
+    def test_T_e_isnan_when_not_set(self):
+        assert np.isnan(self.instance.T_e)
+
+    def test_kappa_isinf_when_not_set(self):
+        assert np.isinf(self.instance.kappa)
diff --git a/plasmapy/atomic/tests/test_ionization_states.py b/plasmapy/atomic/tests/test_ionization_states.py
new file mode 100644
index 00000000..53a81e8e
--- /dev/null
+++ b/plasmapy/atomic/tests/test_ionization_states.py
@@ -0,0 +1,837 @@
+import pytest
+import collections
+from typing import Dict
+from numbers import Real
+import itertools
+
+import numpy as np
+import astropy.units as u
+
+from plasmapy.utils import AtomicError, InvalidIsotopeError, run_test
+from plasmapy.atomic import (
+    State,
+    IonizationState,
+    IonizationStates,
+    atomic_number,
+    mass_number,
+    particle_symbol,
+    Particle,
+)
+
+
+def check_abundances_consistency(abundances: Dict[str, Real], log_abundances: Dict[str, Real]):
+    """
+    Test that a set of abundances is consistent with a set of the base
+    10 logarithm of abundances.
+    """
+    assert abundances.keys() == log_abundances.keys(), (
+        f"Mismatch between keys from abundances and log_abundances.\n\n" 
+        f"abundances.keys():     {abundances.keys()}\n\n"
+        f"log_abundances.keys(): {log_abundances.keys()}")
+
+    for element in abundances.keys():
+        abundance_from_abundances = abundances[element]
+        abundance_from_log_abundances = 10 ** log_abundances[element]
+        assert np.isclose(abundance_from_abundances, abundance_from_log_abundances), (
+            f"Mismatch between abundances and log_abundances.")
+
+def has_attribute(attribute, tests_dict):
+    cases = [test for test in tests_dict.keys() if attribute in tests_dict[test].keys()]
+    if cases:
+        return cases
+    else:
+        raise ValueError(f"No cases with attribute {attribute}")
+
+
+tests = {
+
+    'basic': {
+        'inputs': {'H': [0.1, 0.9], 'helium': (0.2, 0.3, 0.5)},
+        'T_e': 1e6 * u.K,
+        'tol': 1e-15,
+        'abundances': {'hydrogen': 1.0, 'He': 0.1},
+        'kappa': 1.5001,
+    },
+
+    'quantities': {
+        'inputs': {'H': np.array([10, 90]) * u.m ** -3, 'He': np.array([1, 9, 0]) * u.m ** -3},
+    },
+
+    'just H': {
+        'inputs': {'H': [0.1, 0.9]},
+    },
+
+    'H acceptable error': {
+        'inputs': {'H': [1.0, 1e-6]},
+        'tol': 1e-5,
+    },
+
+    'n': {
+        'inputs': {'He': [1, 0, 0], 'H': [1, 0]},
+        'abundances': {'H': 1, 'He': 0.1},
+        'n': 1e9 * u.cm ** -3,
+    },
+
+    'T_e and n': {
+        'inputs': {'H': [0.9, 0.1], 'helium': [0.5, 0.3, 0.2]},
+        'abundances': {'H': 1, 'He': 0.1},
+        'T_e': 1e4 * u.K,
+        'n': 1e15 * u.m ** -3,
+        'kappa': np.inf,
+    },
+
+    'log_abundances': {
+        'inputs': {'H': [1, 0], 'He': [1, 0, 0]},
+        'log_abundances': {'H': 1, 'He': 0},
+        'n': 1e9 * u.cm ** -3,
+    },
+
+    'elements & isotopes': {
+        'inputs': {'H': [0.9, 0.1], 'He-3': [0.3, 0.7, 0.0], 'He-4': [0.29, 0.69, 0.02]},
+        'abundances': {'H': 1, 'He-3': 1e-7, 'He-4': 0.1},
+        'n': 1e12 * u.m ** -3,
+    },
+
+    'ordered elements -> inputs': {
+        'inputs': ['O', 'C', 'H', 'Fe', 'Ar'],
+    },
+
+    'mixed and unordered elements and isotopes': {
+        'inputs': ('Og', 'O', 'H', 'Fe-56', 'He', 'Li-7', 'Li-6'),
+    },
+
+    'number densities -> inputs': {
+        'inputs': {'H': np.array([2, 3]) * u.m ** -3, 'He': np.array([5, 7, 11]) * u.m ** -3},
+    },
+
+    'number densities and n are both inputs': {
+        'inputs': {'H': [0.1, 0.3] * u.cm ** -3},
+        'n': 1e-5 * u.mm ** -3,
+    }
+
+}
+
+test_names = tests.keys()
+
+
+class TestIonizationStates:
+
+    @classmethod
+    def setup_class(cls):
+        cls.instances = {}
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_instantiation(self, test_name):
+        try:
+            self.instances[test_name] = IonizationStates(**tests[test_name])
+        except Exception:
+            pytest.fail(f"Cannot create IonizationStates instance for test='{test_name}'")
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_no_exceptions_from_str(self, test_name):
+        self.instances[test_name].__str__()
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_no_exceptions_from_repr(self, test_name):
+        self.instances[test_name].__repr__()
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_no_exceptions_from_info(self, test_name):
+        self.instances[test_name].info()
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_simple_equality(self, test_name):
+        """Test that __eq__ is not extremely broken."""
+        a = IonizationStates(**tests[test_name])
+        b = IonizationStates(**tests[test_name])
+        assert a == a, f"IonizationStates instance does not equal itself."
+        assert a == b, f"IonizationStates instance does not equal identical instance."
+
+    @pytest.mark.parametrize(
+        'test_name',
+        [test_name for test_name in test_names if isinstance(tests[test_name]['inputs'], dict)],
+    )
+    def test_that_particles_were_set_correctly(self, test_name):
+        input_particles = tests[test_name]['inputs'].keys()
+        particles = [Particle(input_particle) for input_particle in input_particles]
+        expected_particles = {p.particle for p in particles}
+        actual_particles = {particle for particle in self.instances[test_name].ionic_fractions.keys()}
+
+        assert actual_particles == expected_particles, (
+            f"For test='{test_name}', the following should be equal:\n"
+            f"  actual_particles = {actual_particles}\n"
+            f"expected_particles = {expected_particles}")
+
+    @pytest.mark.parametrize('test_name', has_attribute('abundances', tests))
+    def test_that_abundances_kwarg_sets_abundances(self, test_name):
+        try:
+            actual_abundances = self.instances[test_name].abundances
+        except Exception as exc:
+            pytest.fail("Unable to access abundances.")
+
+        elements = set(self.instances[test_name].base_particles)
+        elements_from_abundances = set(actual_abundances.keys())
+
+        if not elements.issubset(elements_from_abundances):
+            pytest.fail(
+                f"The elements whose IonizationStates are being kept "
+                f"track of ({elements}) are not a subset of the "
+                f"elements whose abundances are being kept track of "
+                f"({elements_from_abundances}) for test {test_name}.")
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_that_elements_and_isotopes_are_sorted(self, test_name):
+        elements = self.instances[test_name].base_particles
+        before_sorting = []
+        for element in elements:
+            atomic_numb = atomic_number(element)
+            try:
+                mass_numb = mass_number(element)
+            except InvalidIsotopeError:
+                mass_numb = 0
+            before_sorting.append((atomic_numb, mass_numb))
+        after_sorting = sorted(before_sorting)
+
+        assert before_sorting == after_sorting, (
+            f"Elements/isotopes are not sorted for test='{test_name}':\n"
+            f"  before_sorting = {before_sorting}\n"
+            f"   after_sorting = {after_sorting}\n"
+            f"where above is (atomic_number, mass_number if isotope else 0)")
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_that_ionic_fractions_are_set_correctly(self, test_name):
+
+        errmsg = ""
+
+        elements_actual = self.instances[test_name].base_particles
+        inputs = tests[test_name]["inputs"]
+
+        if isinstance(inputs, dict):
+            input_keys = list(tests[test_name]["inputs"].keys())
+
+            input_keys = sorted(
+                input_keys,
+                key=lambda k: (atomic_number(k), mass_number(k) if Particle(k).isotope else 0)
+            )
+
+            for element, input_key in zip(elements_actual, input_keys):
+                expected = tests[test_name]["inputs"][input_key]
+
+                if isinstance(expected, u.Quantity):
+                    expected = np.array(expected.value / np.sum(expected.value))
+
+                actual = self.instances[test_name].ionic_fractions[element]
+
+                if not np.allclose(actual, expected):
+                    errmsg += (
+                        f"\n\nThere is a discrepancy in ionic fractions for "
+                        f"({test_name}, {element}, {input_key})\n"
+                        f"  expected = {expected}\n"
+                        f"    actual = {actual}")
+
+                if not isinstance(actual, np.ndarray) or isinstance(actual, u.Quantity):
+                    raise AtomicError(f"\n\nNot a numpy.ndarray: ({test_name}, {element})")
+        else:
+            elements_expected = {particle_symbol(element) for element in inputs}
+
+            assert set(self.instances[test_name].base_particles) == elements_expected
+
+            for element in elements_expected:
+                assert all(np.isnan(self.instances[test_name].ionic_fractions[element]))
+        if errmsg:
+            pytest.fail(errmsg)
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_getitem_element(self, test_name):
+        """Test that __get_item__ returns an IonizationState instance"""
+        instance = self.instances[test_name]
+
+        for key in instance.base_particles:
+
+            try:
+                expected = instance.ionic_fractions[key]
+            except Exception as exc:
+                pytest.fail(f"Unable to get ionic_fractions for '{key}' in test='{test_name}'.")
+
+            try:
+                actual = instance[key].ionic_fractions
+            except Exception as exc:
+                pytest(
+                    f"Unable to get item {key} in test={test_name}."
+                )
+
+            try:
+                if all(np.isnan(expected)):
+                    test_passed = True
+                else:
+                    test_passed = np.allclose(expected, actual)
+            except Exception:
+                raise TypeError(
+                    f"For test='{test_name}' and key='{key}', cannot "
+                    f"compare expected ionic fractions of {expected} "
+                    f"with the resulting ionic fractions of {actual}.") from None
+
+            if not test_passed:
+                pytest.fail(
+                    f"For test='{test_name}' and key='{key}', the expected "
+                    f"ionic fractions of {expected} are not all equal "
+                    f"to the resulting ionic fractions of {actual}.")
+
+    @pytest.mark.parametrize('test_name', test_names)
+    def test_getitem_element_intcharge(self, test_name):
+        instance = self.instances[test_name]
+        for particle in instance.base_particles:
+            for int_charge in range(0, atomic_number(particle) + 1):
+                actual = instance[particle, int_charge].ionic_fraction
+                expected = instance.ionic_fractions[particle][int_charge]
+                # We only need to check if one is broken
+            if not np.isnan(actual) and np.isnan(expected):
+                assert np.isclose(actual, expected), (
+                    f"Indexing broken for:\n"
+                    f"       test = '{test_name}'\n"
+                    f"   particle = '{particle}'")
+
+    @pytest.mark.parametrize('test_name', [
+        test_name for test_name in test_names
+        if isinstance(tests[test_name]['inputs'], dict)
+    ])
+    def test_normalization(self, test_name):
+        instance = self.instances[test_name]
+        instance.normalize()
+        not_normalized_elements = []
+        for element in instance.base_particles:
+            is_not_normalized = not np.isclose(
+                np.sum(instance.ionic_fractions[element]),
+                1,
+                atol=1e-19,
+                rtol=0,
+            )
+
+            if is_not_normalized:
+                not_normalized_elements.append(element)
+
+        if not_normalized_elements:
+            pytest.fail(
+                f"In test = '{test_name}', ionic fractions for the "
+                f"following particles were not normalized: "
+                f"{', '.join(not_normalized_elements)}.")
+
+
+def test_abundances_consistency():
+    """Test that ``abundances`` and ``log_abundances`` are consistent."""
+
+    inputs = {'H': [1, 0], 'He': [1, 0, 0]}
+    abundances = {'H': 1.0, 'He': 0.1}
+    elements = abundances.keys()
+
+    log_abundances = {element: np.log10(abundances[element]) for element in elements}
+
+    instance_nolog = IonizationStates(inputs, abundances=abundances)
+    instance_log = IonizationStates(inputs, log_abundances=log_abundances)
+
+    for element in elements:
+        assert np.allclose(
+            instance_log.abundances[element],
+            instance_nolog.abundances[element],
+        ), 'abundances not consistent.'
+
+    for element in elements:
+        assert np.allclose(
+            instance_log.log_abundances[element],
+            instance_nolog.log_abundances[element],
+        ), 'log_abundances not consistent.'
+
+
+class TestIonizationStatesItemAssignment:
+    """
+    Test IonizationStates.__setitem__ and exceptions.
+    """
+
+    @classmethod
+    def setup_class(cls):
+        cls.states = IonizationStates({'H': [0.9, 0.1], 'He': [0.5, 0.4999, 1e-4]})
+
+    @pytest.mark.parametrize("element, new_states", [
+        ('H', [np.nan, np.nan]),
+        ('He', [np.nan, np.nan, np.nan]),
+        ('H', [0.1, 0.9]),
+        ('He', [0.89, 0.1, 0.01]),
+    ])
+    def test_setitem(self, element, new_states):
+        """Test item assignment in an IonizationStates instance."""
+        try:
+            self.states[element] = new_states
+        except Exception:
+            pytest.fail(
+                "Unable to change ionic fractions for an "
+                "IonizationStates instance.")
+        resulting_states = self.states[element].ionic_fractions
+
+        assert np.any([
+            np.allclose(resulting_states, new_states),
+            np.all(np.isnan(resulting_states)) and np.all(np.isnan(new_states)),
+        ])
+
+    @pytest.mark.parametrize('base_particle, new_states, expected_exception', [
+        ('H', (0, 0.9), ValueError),
+        ('H', (-0.1, 1.1), ValueError),
+        ('H', (0.0, 1.0, 0.0), ValueError),
+        ('Li', (0.0, 1.0, 0.0, 0.0), KeyError),
+        ('sdfasd', (0, 1), KeyError),
+        (KeyError, KeyError, KeyError),
+    ])
+    def test_setitem_errors(self, base_particle, new_states, expected_exception):
+        with pytest.raises(expected_exception):
+            self.states[base_particle] = new_states
+
+
+class TestIonizationStatesDensities:
+
+    @classmethod
+    def setup_class(cls):
+
+        cls.initial_ionfracs = {'H': np.array([0.87, 0.13]), 'He': np.array([0.24, 0.37, 0.39])}
+        cls.abundances = {'H': 1.0, 'He': 0.0835}
+        cls.n = 10 * u.m ** -3
+
+        cls.expected_densities = {
+            'H': np.array([8.7, 1.3]) * u.m ** -3,
+            'He': np.array([0.2004, 0.30895, 0.32565]) * u.m ** -3,
+        }
+
+        cls.expected_electron_density = 2.26025 * u.m ** -3
+        cls.states = IonizationStates(cls.initial_ionfracs, abundances=cls.abundances, n=cls.n)
+
+    def test_electron_density(self):
+        assert np.isclose(self.states.n_e.value, self.expected_electron_density.value), (
+            'Mismatch in electron density calculation:\n'
+            f'Calculated = {self.states.n_e}\n'
+            f'Expected   = {self.expected_electron_density}')
+
+    @pytest.mark.parametrize('elem', ['H', 'He'])
+    def test_number_densities(self, elem):
+        assert np.allclose(
+            self.states.number_densities[elem].value,
+            self.expected_densities[elem].value), (
+            f"Mismatch in number densities for {elem}\n"
+            f"Calculated = {self.states.number_densities[elem]}\n"
+            f"Expected   = {self.expected_electron_density}")
+
+
+class TestIonizationStatesAttributes:
+
+    @classmethod
+    def setup_class(cls):
+        cls.elements = ['H', 'He', 'Li', 'Fe']
+        cls.instance = IonizationStates(cls.elements)
+        cls.new_n = 5.153 * u.cm ** -3
+
+    @pytest.mark.parametrize("uninitialized_attribute", ['T_e', 'n', 'n_e'])
+    def test_attribute_defaults_to_nan(self, uninitialized_attribute):
+        command = f"self.instance.{uninitialized_attribute}"
+        default_value = eval(command)
+        assert np.isnan(default_value), (
+            f"{uninitialized_attribute} does not default to nan but "
+            f"instead defaults to {default_value}.")
+
+    def test_kappa_defaults_to_inf(self):
+        assert np.isinf(self.instance.kappa), \
+            "kappa does not default to a value of inf."
+
+    @pytest.mark.parametrize(
+        "uninitialized_attribute",
+        ["number_densities", "ionic_fractions"])
+    def test_attribute_defaults_to_dict_of_nans(self, uninitialized_attribute):
+        command = f"self.instance.{uninitialized_attribute}"
+        default_value = eval(command)
+        assert list(default_value.keys()) == self.elements, "Incorrect base particle keys."
+        for element in self.elements:
+            assert len(default_value[element]) == atomic_number(element) + 1, \
+                f"Incorrect number of ionization levels for {element}."
+            assert np.all(np.isnan(default_value[element])), (
+                f"The values do not default to an array of nans for "
+                f"{element}.")
+
+    @pytest.mark.parametrize("uninitialized_attribute", ["abundances", "log_abundances"])
+    def test_abundances_default_to_nans(self, uninitialized_attribute):
+        command = f"self.instance.{uninitialized_attribute}"
+        default_value = eval(command)
+        for element in self.elements:
+            assert isinstance(default_value[element], Real)
+            assert np.isnan(default_value[element])
+
+    @pytest.mark.parametrize(
+        "attribute, invalid_value, expected_exception", [
+            ('T_e', '5 * u.m', u.UnitsError),
+            ('T_e', '-1 * u.K', AtomicError),
+            ('n', '5 * u.m', u.UnitsError),
+            ('n', '-1 * u.m ** -3', AtomicError),
+            ('ionic_fractions', {'H': [0.3, 0.7], 'He': [-0.1, 0.4, 0.7]}, AtomicError),
+            ('ionic_fractions', {'H': [0.3, 0.7], 'He': [1.01, 0.0, 0.7]}, AtomicError),
+            ('ionic_fractions', {'H': [0.3, 0.6], 'He': [1.0, 0.0, 0.0]}, AtomicError),
+            ('ionic_fractions', {'H': [1.0, 0.0]}, AtomicError),
+        ])
+    def test_attribute_exceptions(self, attribute, invalid_value, expected_exception):
+
+        command = f"self.instance.{attribute} = {invalid_value}"
+        errmsg = f"No {expected_exception} was raised for command\n\n: {command}"
+
+        with pytest.raises(expected_exception, message=errmsg):
+            exec(command)
+
+    def test_setting_ionic_fractions_for_single_element(self):
+        """
+        Test that __setitem__ correctly sets new ionic fractions when
+        used for just H, while not changing the ionic fractions for He
+        from the uninitialized default of an array of nans of length 3.
+        """
+        self.new_fractions = [0.3, 0.7]
+        self.instance['H'] = self.new_fractions
+        resulting_fractions = self.instance.ionic_fractions['H']
+        assert np.allclose(self.new_fractions, resulting_fractions), \
+            "Ionic fractions for H not set using __setitem__."
+        assert 'He' in self.instance.ionic_fractions.keys(), (
+            "He is missing in ionic_fractions after __setitem__ was "
+            "used to set H ionic fractions.")
+        assert np.all(np.isnan(self.instance.ionic_fractions['He'])), (
+            "He ionic fractions are not all nans after __setitem__ "
+            "was used to set H ionic fractions.")
+
+    @pytest.mark.parametrize(
+        "key, invalid_fracs, expected_exception",[
+            ('H', [-0.01, 1.01], ValueError),
+            ('H', [0.4, 0.5], ValueError),
+            ('H', [0.5, 0.5, 0.0], ValueError),
+            ('He', [0.5, 0.5], ValueError),
+            ('He', [0.1, 0.9, 0.0, 0.0], ValueError),
+            ('He', [0.9, 0.1, 0.0] * u.m ** -2, ValueError),
+            ('He', [-0.01, 0.99, 0.02], ValueError),
+            ('He', [1.01, -0.02, 0.01], ValueError),
+        ])
+    def test_setting_invalid_ionfracs(self, key, invalid_fracs, expected_exception):
+        errmsg = (
+            f"No {expected_exception} is raised when trying to assign "
+            f"{invalid_fracs} to {key} in an IonizationStates instance.")
+        with pytest.raises(expected_exception, message=errmsg):
+            self.instance[key] = invalid_fracs
+
+    def test_setting_incomplete_abundances(self):
+        new_abundances = {'H': 1, 'He': 0.1, 'Fe': 1e-5, 'Au': 1e-8}  # missing lithium
+        with pytest.raises(AtomicError):
+            self.instance.abundances = new_abundances
+
+    def test_setting_abundances(self):
+        new_abundances = {'H': 1, 'He': 0.1, 'Li': 1e-4, 'Fe': 1e-5, 'Au': 1e-8}
+
+        log_new_abundances = {
+            element: np.log10(new_abundances[element])
+            for element in new_abundances.keys()
+        }
+
+        try:
+            self.instance.abundances = new_abundances
+        except Exception:
+            pytest(f"Could not set abundances to {new_abundances}.")
+        else:
+            check_abundances_consistency(self.instance.abundances, self.instance.log_abundances)
+
+        try:
+            self.instance.log_abundances = log_new_abundances
+        except Exception:
+            pytest.fail(f"Could not set log_abundances to {log_new_abundances}.")
+        else:
+            check_abundances_consistency(self.instance.abundances, self.instance.log_abundances)
+
+    @pytest.mark.parametrize(
+        "invalid_indices", [
+            (1, 2, 3),
+            'C',
+            'H-1',
+            ('Fe', -1),
+            ('Fe', 27),
+            ('He', -1),
+            ('He', 3),
+            ('Fe', slice(3, 7)),
+        ])
+    def test_invalid_indices(self, invalid_indices):
+        with pytest.raises(IndexError):
+            self.instance[invalid_indices]
+
+    @pytest.mark.parametrize("index", ["H", "Fe"])
+    def test_getitem_one_index(self, index):
+        instance = self.instance
+        result = instance[index]
+
+        if np.all(np.isnan(instance.number_densities[index])):
+            inputs = instance.ionic_fractions[index]
+        else:
+            inputs = instance.number_densities[index]
+
+        expected = IonizationState(index, inputs, T_e=instance.T_e, kappa=instance.kappa)
+
+        assert isinstance(result, IonizationState)
+        assert result == expected
+
+    @pytest.mark.parametrize("indices", [("H", 1), ("Fe", 6)])
+    def test_getitem_two_indices(self, indices):
+        instance = self.instance
+        result = instance[indices]
+
+        particle = indices[0]
+        integer_charge = indices[1]
+
+        assert isinstance(result, State)
+        assert result.integer_charge == integer_charge
+
+        expected_ionic_fraction = instance.ionic_fractions[particle][integer_charge]
+
+        assert np.any([
+            np.isclose(result.ionic_fraction, expected_ionic_fraction),
+            np.isnan(result.ionic_fraction) and np.isnan(expected_ionic_fraction),
+        ])
+
+        assert result.ionic_symbol == particle_symbol(particle, Z=integer_charge)
+
+    def test_setting_n(self):
+        try:
+            self.instance.n = self.new_n
+        except Exception:
+            pytest.fail("Unable to set number density scaling factor attribute")
+        if not u.quantity.allclose(self.instance.n, self.new_n):
+            pytest.fail("Number density scaling factor was not set correctly.")
+        if not self.instance.n.unit == u.m ** -3:
+            pytest.fail("Incorrect units for new number density.")
+
+    def test_resetting_valid_densities(self):
+        """
+        Test that item assignment can be used to set number densities
+        that preserve the total element number density.
+        """
+
+        element = "H"
+        valid_ionic_fractions = [0.54, 0.46]
+        original_n_elem = np.sum(self.instance.number_densities[element])
+        valid_number_densities = valid_ionic_fractions * original_n_elem
+
+        try:
+            self.instance[element] = valid_number_densities
+        except Exception:
+            pytest.fail("Unable to set valid number densities using item assignment.")
+
+        assert u.quantity.allclose(
+            self.instance.ionic_fractions[element],
+            valid_ionic_fractions,
+        ), "Item assignment of valid number densities did not yield correct ionic fractions."
+
+        assert u.quantity.allclose(
+            self.instance.number_densities[element],
+            valid_number_densities,
+        ), "Item assignment of valid number densities did not yield correct number densities."
+
+    def test_resetting_invalid_densities(self):
+        """
+        Test that item assignment with number densities that would
+        change the total element number density raises an exception.
+        """
+        element = "H"
+        original_n_elem = np.sum(self.instance.number_densities[element])
+        invalid_number_densities = np.array([1.0001, 0]) * original_n_elem
+        with pytest.raises(ValueError):
+            self.instance[element] = invalid_number_densities
+
+    def test_elemental_abundances_not_quantities(self):
+        for element in self.instance.base_particles:
+            assert not isinstance(self.instance.abundances[element], u.Quantity)
+
+    @pytest.mark.parametrize("element", ["H", "He", "Fe"])
+    def test_ionic_fractions_not_quantities(self, element):
+        ionic_fractions = self.instance.ionic_fractions[element]
+        if isinstance(ionic_fractions, u.Quantity):
+            pytest.fail(f"The ionic fractions of {element} are a Quantity but should not be.")
+
+    def test_that_iron_ionic_fractions_are_still_undefined(self):
+        assert 'Fe' in self.instance.ionic_fractions.keys()
+        iron_fractions = self.instance.ionic_fractions['Fe']
+        assert len(iron_fractions) == atomic_number('Fe') + 1
+        assert np.all(np.isnan(iron_fractions))
+
+    def test_base_particles(self):
+        """
+        Test that the original base particles remain as the base
+        particles after performing a bunch of operations that should not
+        change them.
+        """
+        assert self.instance.base_particles == self.elements
+
+    def test_base_particles_equal_ionic_fraction_particles(self):
+        assert self.instance.base_particles == list(self.instance.ionic_fractions.keys())
+
+
+IE = collections.namedtuple("IE", ["inputs", "expected_exception"])
+
+tests_for_exceptions = {
+
+    'wrong type': IE({
+        "inputs": None,
+    }, AtomicError),
+
+    'not normalized': IE({
+        "inputs": {'He': [0.4, 0.5, 0.0]}, "tol": 1e-9
+    }, AtomicError),
+
+    'negative ionfrac': IE({
+        "inputs": {'H': [-0.1, 1.1]},
+    }, AtomicError),
+
+    'ion': IE({
+        "inputs": {'H': [0.1, 0.9], 'He+': [0.0, 0.9, 0.1]},
+    }, AtomicError),
+
+    'repeat elements': IE({
+        "inputs": {'H': [0.1, 0.9], "hydrogen": [0.2, 0.8]},
+    }, AtomicError),
+
+    'isotope of element': IE({
+        "inputs": {'H': [0.1, 0.9], "D": [0.2, 0.8]},
+    }, AtomicError),
+
+    'negative abundance': IE({
+        "inputs": {"H": [0.1, 0.9], "He": [0.4, 0.5, 0.1]},
+        "abundances": {"H": 1, "He": -0.1},
+    }, AtomicError),
+
+    'imaginary abundance': IE({
+        "inputs": {"H": [0.1, 0.9], "He": [0.4, 0.5, 0.1]},
+        "abundances": {"H": 1, "He": 0.1j},
+    }, AtomicError),
+
+    'wrong density units': IE({
+        "inputs": {"H": [10, 90] * u.m ** -3, "He": [0.1, 0.9, 0] * u.m ** -2},
+        "abundances": {"H": 1, "He": 0.1},
+    }, AtomicError),
+
+    'abundance redundance': IE({
+        "inputs": {"H": [10, 90] * u.m ** -3, "He": [0.1, 0.9, 0] * u.m ** -3},
+        "abundances": {"H": 1, "He": 0.1},
+    }, AtomicError),
+
+    'abundance contradiction': IE({
+        "inputs": {"H": [10, 90] * u.m ** -3, "He": [0.1, 0.9, 0] * u.m ** -3},
+        "abundances": {"H": 1, "He": 0.11},
+    }, AtomicError),
+
+    'kappa too small': IE({
+        "inputs": ["H"],
+        "kappa": 1.499999,
+    }, AtomicError),
+
+    'negative n': IE({
+        "inputs": ["H"],
+        "n": -1 * u.cm ** -3,
+    }, AtomicError),
+
+    'negative T_e': IE({
+        "inputs": ["H-1"],
+        "T_e": -1 * u.K,
+    }, AtomicError),
+}
+
+
+@pytest.mark.parametrize('test_name', tests_for_exceptions.keys())
+def test_exceptions_upon_instantiation(test_name):
+    """
+    Test that appropriate exceptions are raised for inappropriate inputs
+    to IonizationStates when first instantiated.
+    """
+    run_test(
+        IonizationStates,
+        kwargs=tests_for_exceptions[test_name].inputs,
+        expected_outcome=tests_for_exceptions[test_name].expected_exception,
+    )
+
+
+class TestIonizationStatesDensityEqualities:
+    """
+    Test that IonizationStates instances are equal or not equal to each
+    other as they should be for different combinations of inputs
+    related to ionic_fractions, number densities, and abundances.
+    """
+
+    @classmethod
+    def setup_class(cls):
+
+        # Create arguments to IonizationStates that are all consistent
+        # with each other.
+
+        cls.ionic_fractions = {'H': [0.9, 0.1], 'He': [0.99, 0.01, 0.0]}
+        cls.abundances = {'H': 1, 'He': 0.08}
+        cls.n = 5.3 * u.m ** -3
+        cls.number_densities = {
+            element: cls.ionic_fractions[element] * cls.n * cls.abundances[element]
+            for element in cls.ionic_fractions.keys()
+        }
+
+        # The keys that begin with 'ndens' have enough information to
+        # yield the number_densities attribute, whereas the keys that
+        # begin with "no_ndens" do not.
+
+        cls.dict_of_kwargs = {
+            'ndens1': {'inputs': cls.ionic_fractions, 'abundances': cls.abundances, 'n': cls.n},
+            'ndens2': {'inputs': cls.number_densities},
+            'no_ndens3': {'inputs': cls.ionic_fractions},
+            'no_ndens4': {'inputs': cls.ionic_fractions, 'abundances': cls.abundances},
+            'no_ndens5': {'inputs': cls.ionic_fractions, 'n': cls.n},
+        }
+
+        cls.instances = {
+            key: IonizationStates(**cls.dict_of_kwargs[key])
+            for key in cls.dict_of_kwargs.keys()
+        }
+
+    @pytest.mark.parametrize("test_key", ['ndens1', 'ndens2'])
+    def test_number_densities_defined(self, test_key):
+        number_densities = self.instances[test_key].number_densities
+        for base_particle in self.instances[test_key].base_particles:
+            assert not np.any(np.isnan(number_densities[base_particle])), (
+                f"Test {test_key} should have number densities "
+                f"defined, but doesn't.")
+
+    @pytest.mark.parametrize("test_key", ['no_ndens3', 'no_ndens4', 'no_ndens5'])
+    def test_number_densities_undefined(self, test_key):
+        number_densities = self.instances[test_key].number_densities
+        for base_particle in self.instances[test_key].base_particles:
+            assert np.all(np.isnan(number_densities[base_particle])), (
+                f"Test {test_key} should not have number densities "
+                f"defined, but does.")
+
+    @pytest.mark.parametrize("this, that", itertools.product(
+        ['ndens1', 'ndens2', 'no_ndens3', 'no_ndens4', 'no_ndens5'],
+        repeat=2,
+    ))
+    def test_equality(self, this, that):
+        """
+        Test that the IonizationStates instances that should provide
+        ``number_densities`` are all equal to each other.  Test that the
+        instances that should not provide ``number_densities`` are all
+        equal to each other.  Test that each instance that should
+        provide ``number_densities`` is not equal to each instance that
+        should not provide ``number_densities``.
+        """
+        expect_equality = this[0:4] == that[0:4]
+        are_equal = self.instances[this] == self.instances[that]
+        if expect_equality != are_equal:
+            print(f"{this} kwargs:\n {self.dict_of_kwargs[this]}\n")
+            self.instances[this].info()
+            print()
+            print(f"{that} kwargs:\n {self.dict_of_kwargs[that]}\n")
+            self.instances[that].info()
+            descriptor = "equal" if expect_equality else "unequal"
+            pytest.fail(
+                f"Cases {this} and {that} should be {descriptor} but "
+                f"are not.")
+
+
+def test_number_density_assignment():
+    instance = IonizationStates(["H", "He"])
+    number_densities = [2, 3, 5] * u.m ** -3
+    instance["He"] = number_densities
diff --git a/plasmapy/atomic/tests/test_particle_class.py b/plasmapy/atomic/tests/test_particle_class.py
index 59b282f1..881e74e0 100644
--- a/plasmapy/atomic/tests/test_particle_class.py
+++ b/plasmapy/atomic/tests/test_particle_class.py
@@ -1,12 +1,17 @@
 import pytest
-import numpy as np
-from astropy import units as u
 import inspect
-from plasmapy.utils import roman
+import collections
 
-from ...constants import m_p, m_e, m_n, e, c
+import numpy as np
+from astropy import units as u
 
-from ...utils import (
+from plasmapy.utils import roman
+from plasmapy.constants import m_p, m_e, m_n, e, c
+from plasmapy.atomic.atomic import known_isotopes
+from plasmapy.atomic.isotopes import _Isotopes
+from plasmapy.atomic.particle_class import Particle
+from plasmapy.atomic.special_particles import ParticleZoo
+from plasmapy.utils import (
     AtomicWarning,
     AtomicError,
     MissingAtomicDataError,
@@ -20,12 +25,7 @@
     run_test_equivalent_calls,
 )
 
-from ..atomic import known_isotopes
-from ..isotopes import _Isotopes
-from ..particle_class import Particle
-from ..special_particles import ParticleZoo
-
-# (arg, kwargs, results_dict
+# (arg, kwargs, results_dict)
 test_Particle_table = [
 
     ('neutron', {},
@@ -559,7 +559,7 @@ def test_Particle_warnings(arg, kwargs, attribute, warning):
 
 
 def test_Particle_cmp():
-    """Test __eq__ and __ne__ in the Particle class."""
+    """Test ``__eq__`` and ``__ne__`` in the Particle class."""
     proton1 = Particle('p+')
     proton2 = Particle('proton')
     electron = Particle('e-')
@@ -589,7 +589,7 @@ def test_Particle_cmp():
 @pytest.mark.parametrize('isotope, ion', nuclide_mass_and_mass_equiv_table)
 def test_particle_class_mass_nuclide_mass(isotope: str, ion: str):
     """
-    Test that the `mass` and `nuclide_mass` attributes return
+    Test that the ``mass`` and ``nuclide_mass`` attributes return
     equivalent values when appropriate.  The inputs should generally be
     an isotope with no charge information, and a fully ionized ion of
     that isotope, in order to make sure that the nuclide mass of the
@@ -610,12 +610,13 @@ def test_particle_class_mass_nuclide_mass(isotope: str, ion: str):
 
     else:
 
-        inputerrmsg = (f"isotope = {repr(isotope)} and ion = {repr(ion)} are "
-                       f"not valid inputs to this test. The inputs should be "
-                       f"an isotope with no charge information, and a fully "
-                       f"ionized ion of that isotope, in order to make sure "
-                       f"that the nuclide mass of the isotope equals the mass "
-                       f"of the ion.")
+        inputerrmsg = (
+            f"isotope = {repr(isotope)} and ion = {repr(ion)} are "
+            f"not valid inputs to this test. The inputs should be "
+            f"an isotope with no charge information, and a fully "
+            f"ionized ion of that isotope, in order to make sure "
+            f"that the nuclide mass of the isotope equals the mass "
+            f"of the ion.")
 
         assert Isotope.isotope and not Isotope.ion, inputerrmsg
         assert Isotope.isotope == Ion.isotope, inputerrmsg
@@ -634,7 +635,7 @@ def test_particle_half_life_string():
     """
     Find the first isotope where the half-life is stored as a string
     (because the uncertainties are too great), and tests that requesting
-    the half-life of that isotope causes a MissingAtomicDataWarning
+    the half-life of that isotope causes a `MissingAtomicDataWarning`
     whilst returning a string.
     """
 
@@ -701,8 +702,8 @@ def opposite(particle):
         opposite_particle = ~particle
     except Exception as exc:
         raise InvalidParticleError(
-                f"The unary ~ (invert) operator is unable to find the "
-                f"antiparticle of {particle}.") from exc
+            f"The unary ~ (invert) operator is unable to find the "
+            f"antiparticle of {particle}.") from exc
     return opposite_particle
 
 
@@ -769,3 +770,41 @@ def test_particleing_a_particle(arg):
         f"Particle({repr(arg)}) is the same object in memory as "
         f"Particle(Particle({repr(arg)})), when it is intended to "
         f"create a new object in memory (e.g., a copy).")
+
+
+@pytest.mark.parametrize("key", [Particle('H'), Particle('e+')])
+def test_that_object_can_be_dict_key(key):
+    """
+    Test that ``key`` can be used as the key of a `dict`.
+
+    This test will fail if ``key`` does not equal itself or if ``key``
+    is not hashable.  If this test fails, there is a problem with the
+    ``__eq__`` and ``__hash__`` methods of ``key``.
+
+    In most cases, objects that are mutable should not be hashable since
+    they may change.
+
+    """
+    # TODO: I wrote this to be pretty general since I felt like
+    # TODO: procrastinating other things, so we can probably put this
+    # TODO: into utils.pytest_helpers later on.  There are likely other
+    # TODO: classes that should be able to be used as keys of dicts.
+
+    value = 42
+
+    try:
+        dictionary = {key: value}
+    except Exception as exc:
+        error_message = f"{key} is not a valid key for a dict. "
+        if not isinstance(key, collections.Hashable):
+            error_message += f"{key} is not hashable. "
+        try:
+            key_equals_itself = key == key
+        except Exception:
+            error_message += f"{key} == {key} cannot be evaluated. "
+        else:
+            if not key_equals_itself:
+                error_message += f"{key} does not equal itself."
+        raise TypeError(error_message) from exc
+
+    assert dictionary[key] is value
diff --git a/plasmapy/physics/tests/test_parameters.py b/plasmapy/physics/tests/test_parameters.py
index 63686127..03ca5588 100644
--- a/plasmapy/physics/tests/test_parameters.py
+++ b/plasmapy/physics/tests/test_parameters.py
@@ -809,28 +809,28 @@ def test_inertial_length():
     with pytest.raises(ValueError):
         inertial_length(-5 * u.m ** -3, particle='p')
 
-    with pytest.raises(ValueError):
+    with pytest.raises(InvalidParticleError):
         inertial_length(n_i, particle=-135)
 
     with pytest.warns(u.UnitsWarning):
         inertial_length_no_units = inertial_length(1e19, particle='p')
         assert inertial_length_no_units == inertial_length(1e19 * u.m ** -3, particle='p')
 
-    assert inertial_length(n_e).unit.is_equivalent(u.m)
+    assert inertial_length(n_e, 'e-').unit.is_equivalent(u.m)
 
-    assert np.isclose(inertial_length(1 * u.cm ** -3).cgs.value, 5.31e5, rtol=1e-3)
+    assert np.isclose(inertial_length(1 * u.cm ** -3, 'e-').cgs.value, 5.31e5, rtol=1e-3)
 
     with pytest.warns(u.UnitsWarning):
-        inertial_length(5)
+        inertial_length(5, 'e-')
 
     with pytest.raises(u.UnitConversionError):
-        inertial_length(5 * u.m)
+        inertial_length(5 * u.m, 'e-')
 
     with pytest.raises(ValueError):
-        inertial_length(-5 * u.m ** -3)
+        inertial_length(-5 * u.m ** -3, 'e-')
 
     with pytest.warns(u.UnitsWarning):
-        assert inertial_length(1e19) == inertial_length(1e19 * u.m ** -3)
+        assert inertial_length(1e19, 'e-') == inertial_length(1e19 * u.m ** -3, 'e-')
 
     assert_can_handle_nparray(inertial_length)