Repository: PlasmaPy/PlasmaPy
Pull request #404 (2018-04-24 22:31:49)

Create classes to represent ionization state distributions
My plan for this PR is to create classes to represent the ionization state distributions of one or more elements.  I am going to add in a bunch of dunder methods like `__getitem__` and maybe `__call__` to help making access to the ionization states more straightfoward and intuitive.  Any suggestions on the naming convention will be helpful so that we can maximize readability. 

Eventually we'll need a way to calculate ionization state distributions assuming collisional ionization equilibrium, but that will be for a different PR.  The purpose of this PR is to set up how to store and access the ionization distributions.  This will be discussed in #352.

This will address some of #352.  It will probably be best to wait until after the `0.1.0` release to merge this, since this PR is only for a partial implementation anyway.



pep8speaks: Hello @namurphy! Thanks for submitting your pull request.



 - In the file [`plasmapy/atomic/ionization_states.py`](https://github.com/PlasmaPy/PlasmaPy/blob/a37e283c255b8caea6b96e5fac6368dc122a6bca/plasmapy/atomic/ionization_states.py), following are the PEP8 issues :

> [Line 6:1](https://github.com/PlasmaPy/PlasmaPy/blob/a37e283c255b8caea6b96e5fac6368dc122a6bca/plasmapy/atomic/ionization_states.py#L6): [E302](https://duckduckgo.com/?q=pep8%20E302) expected 2 blank lines, found 1
> [Line 65:5](https://github.com/PlasmaPy/PlasmaPy/blob/a37e283c255b8caea6b96e5fac6368dc122a6bca/plasmapy/atomic/ionization_states.py#L65): [E303](https://duckduckgo.com/?q=pep8%20E303) too many blank lines (3)
> [Line 87:22](https://github.com/PlasmaPy/PlasmaPy/blob/a37e283c255b8caea6b96e5fac6368dc122a6bca/plasmapy/atomic/ionization_states.py#L87): [W292](https://duckduckgo.com/?q=pep8%20W292) no newline at end of file



PlasmaPy follows the [PEP8 style guide](https://www.python.org/dev/peps/pep-0008/), which may be checked locally using [pycodestyle](http://pycodestyle.pycqa.org/en/latest/intro.html)
codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=h1) Report
> Merging [#404](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/87515b4a27052001e11213cc97341dbd9a2f3590?src=pr&el=desc) will **decrease** coverage by `1.39%`.
> The diff coverage is `78.89%`.

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

```diff
@@            Coverage Diff            @@
##           master     #404     +/-   ##
=========================================
- Coverage   92.85%   91.45%   -1.4%     
=========================================
  Files          67       70      +3     
  Lines        6434     7128    +694     
=========================================
+ Hits         5974     6519    +545     
- Misses        460      609    +149
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/atomic/particle\_class.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL3BhcnRpY2xlX2NsYXNzLnB5) | `97.18% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/parsing.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL3BhcnNpbmcucHk=) | `99.46% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/atomic.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2F0b21pYy5weQ==) | `93.33% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL19faW5pdF9fLnB5) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/particle\_input.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL3BhcnRpY2xlX2lucHV0LnB5) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/ionization\_states.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGVzLnB5) | `74.04% <74.04%> (ø)` | |
| [plasmapy/atomic/tests/test\_ionization\_states.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL3Rlc3RzL3Rlc3RfaW9uaXphdGlvbl9zdGF0ZXMucHk=) | `85.71% <85.71%> (ø)` | |
| [plasmapy/atomic/tests/test\_ionization\_state.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL3Rlc3RzL3Rlc3RfaW9uaXphdGlvbl9zdGF0ZS5weQ==) | `87.5% <87.5%> (ø)` | |
| ... and [1 more](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree-more) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?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/404?src=pr&el=footer). Last update [87515b4...89f8861](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

namurphy: Thanks for the code review, @StanczakDominik!  We've been using this and related functionality in a separate but related [project for non-equilibrium ionization (NEI) modeling](https://github.com/NEI-modeling/NEI) that we're still developing and refining, and that I'm hoping will get implemented into a different project.  

There are a few things I've been wondering about on larger scales.  In which project should the `IonizationStates` class (or something similar) go?  And also...where should NEI modeling functionality go?  The ionization state of a plasma is a really basic property, so I've been thinking that it should go in PlasmaPy, even though fiasco and SpectroscoPyx are also really good possibilities.  NEI modeling is a bit more specialized, but the data on recombination and ionization rates is pretty much the same data that's needed to calculate equilibrium ionization fractions, which is a really fundamental thing too so PlasmaPy might be the best place.  

I should probably create a separate issue to ask this question, but I'm going on vacation in three minutes so it'll have to wait! :bullettrain_front: :bridge_at_night: 

codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=h1) Report
> Merging [#404](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/0c6fa440cb8c8fea4c5ee70142c009525585e221?src=pr&el=desc) will **decrease** coverage by `2.2%`.
> The diff coverage is `79.03%`.

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

```diff
@@            Coverage Diff             @@
##           master     #404      +/-   ##
==========================================
- Coverage   98.07%   95.87%   -2.21%     
==========================================
  Files          46       48       +2     
  Lines        3796     4291     +495     
==========================================
+ Hits         3723     4114     +391     
- Misses         73      177     +104
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/atomic/particle\_class.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL3BhcnRpY2xlX2NsYXNzLnB5) | `100% <ø> (ø)` | :arrow_up: |
| [plasmapy/atomic/atomic.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2F0b21pYy5weQ==) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL19faW5pdF9fLnB5) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/ionization\_states.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGVzLnB5) | `69.89% <69.89%> (ø)` | |
| [plasmapy/atomic/ionization\_state.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGUucHk=) | `90.69% <90.69%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?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/404?src=pr&el=footer). Last update [0c6fa44...22dc045](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=h1) Report
> Merging [#404](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/0c6fa440cb8c8fea4c5ee70142c009525585e221?src=pr&el=desc) will **decrease** coverage by `2.2%`.
> The diff coverage is `79.03%`.

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

```diff
@@            Coverage Diff             @@
##           master     #404      +/-   ##
==========================================
- Coverage   98.07%   95.87%   -2.21%     
==========================================
  Files          46       48       +2     
  Lines        3796     4291     +495     
==========================================
+ Hits         3723     4114     +391     
- Misses         73      177     +104
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/atomic/particle\_class.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL3BhcnRpY2xlX2NsYXNzLnB5) | `100% <ø> (ø)` | :arrow_up: |
| [plasmapy/atomic/atomic.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2F0b21pYy5weQ==) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL19faW5pdF9fLnB5) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/ionization\_states.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGVzLnB5) | `69.89% <69.89%> (ø)` | |
| [plasmapy/atomic/ionization\_state.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGUucHk=) | `90.69% <90.69%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?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/404?src=pr&el=footer). Last update [0c6fa44...913b204](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=h1) Report
> Merging [#404](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/0c6fa440cb8c8fea4c5ee70142c009525585e221?src=pr&el=desc) will **decrease** coverage by `2.35%`.
> The diff coverage is `78.23%`.

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

```diff
@@            Coverage Diff             @@
##           master     #404      +/-   ##
==========================================
- Coverage   98.07%   95.72%   -2.36%     
==========================================
  Files          46       48       +2     
  Lines        3796     4305     +509     
==========================================
+ Hits         3723     4121     +398     
- Misses         73      184     +111
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/atomic/particle\_class.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL3BhcnRpY2xlX2NsYXNzLnB5) | `100% <ø> (ø)` | :arrow_up: |
| [plasmapy/atomic/atomic.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2F0b21pYy5weQ==) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL19faW5pdF9fLnB5) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/ionization\_states.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGVzLnB5) | `69.89% <69.89%> (ø)` | |
| [plasmapy/atomic/ionization\_state.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGUucHk=) | `88.2% <88.2%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?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/404?src=pr&el=footer). Last update [0c6fa44...913b204](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=h1) Report
> Merging [#404](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/0c6fa440cb8c8fea4c5ee70142c009525585e221?src=pr&el=desc) will **decrease** coverage by `2.35%`.
> The diff coverage is `78.23%`.

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

```diff
@@            Coverage Diff             @@
##           master     #404      +/-   ##
==========================================
- Coverage   98.07%   95.72%   -2.36%     
==========================================
  Files          46       48       +2     
  Lines        3796     4305     +509     
==========================================
+ Hits         3723     4121     +398     
- Misses         73      184     +111
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/atomic/particle\_class.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL3BhcnRpY2xlX2NsYXNzLnB5) | `100% <ø> (ø)` | :arrow_up: |
| [plasmapy/atomic/atomic.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2F0b21pYy5weQ==) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL19faW5pdF9fLnB5) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/ionization\_states.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGVzLnB5) | `69.89% <69.89%> (ø)` | |
| [plasmapy/atomic/ionization\_state.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGUucHk=) | `88.2% <88.2%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?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/404?src=pr&el=footer). Last update [0c6fa44...913b204](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=h1) Report
> Merging [#404](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/0c6fa440cb8c8fea4c5ee70142c009525585e221?src=pr&el=desc) will **decrease** coverage by `2.35%`.
> The diff coverage is `78.23%`.

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

```diff
@@            Coverage Diff             @@
##           master     #404      +/-   ##
==========================================
- Coverage   98.07%   95.72%   -2.36%     
==========================================
  Files          46       48       +2     
  Lines        3796     4305     +509     
==========================================
+ Hits         3723     4121     +398     
- Misses         73      184     +111
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/atomic/particle\_class.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL3BhcnRpY2xlX2NsYXNzLnB5) | `100% <ø> (ø)` | :arrow_up: |
| [plasmapy/atomic/atomic.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2F0b21pYy5weQ==) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL19faW5pdF9fLnB5) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/ionization\_states.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGVzLnB5) | `69.89% <69.89%> (ø)` | |
| [plasmapy/atomic/ionization\_state.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGUucHk=) | `88.2% <88.2%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?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/404?src=pr&el=footer). Last update [0c6fa44...913b204](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=h1) Report
> Merging [#404](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/0c6fa440cb8c8fea4c5ee70142c009525585e221?src=pr&el=desc) will **decrease** coverage by `2.35%`.
> The diff coverage is `78.23%`.

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

```diff
@@            Coverage Diff             @@
##           master     #404      +/-   ##
==========================================
- Coverage   98.07%   95.72%   -2.36%     
==========================================
  Files          46       48       +2     
  Lines        3796     4305     +509     
==========================================
+ Hits         3723     4121     +398     
- Misses         73      184     +111
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/atomic/particle\_class.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL3BhcnRpY2xlX2NsYXNzLnB5) | `100% <ø> (ø)` | :arrow_up: |
| [plasmapy/atomic/atomic.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2F0b21pYy5weQ==) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL19faW5pdF9fLnB5) | `100% <100%> (ø)` | :arrow_up: |
| [plasmapy/atomic/ionization\_states.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGVzLnB5) | `69.89% <69.89%> (ø)` | |
| [plasmapy/atomic/ionization\_state.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404/diff?src=pr&el=tree#diff-cGxhc21hcHkvYXRvbWljL2lvbml6YXRpb25fc3RhdGUucHk=) | `88.2% <88.2%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?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/404?src=pr&el=footer). Last update [0c6fa44...25d073a](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/404?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

StanczakDominik: Also, to make sure this doesn't get lost in GitHub's review interface, [this comment discusses the test failure](https://github.com/PlasmaPy/PlasmaPy/pull/404#discussion_r224181200).
namurphy: > **Transport PR Incident**

:joy_cat: 

I still have more tests and docstrings to put in, but I don't have to finish all of them in this PR.  I was planning on working on this for one more day when I can take some time away from Fortran, and then start working on #521.  This is the last of my backlog of PRs to finish up, so I'm eager to get it done soon!
StanczakDominik: The Travis builds are only failing because I messed up a few git remote magic tricks; they're passing [over here](https://travis-ci.org/PlasmaPy/PlasmaPy/builds/441512524).
namurphy: The downside of having `__repr__` return multiple lines of diagnostic data is that we can end up with messy error messages with `'\n'` showing up a bunch of times.
```
___________ Test_IonizationStates.test_inequality[number densities -> inputs] ___________
plasmapy/atomic/tests/test_ionization_states.py:159: in test_inequality
    assert original_instance != altered_instance
E   assert IonizationStates instance for: H, He\n\n  H  0+: 0.400\n  H  1+: 0.600\n\n  He  0+: 0.217\n  He  1+: 0.304\n  He  2+: 0.478\n\n != IonizationStates instance for: H, He\n\n  H  0+: 0.400\n  H  1+: 0.600\n\n  He  0+: 0.217\n  He  1+: 0.304\n  He  2+: 0.478\n\n
```
I'm still wondering if being able to easily print out really useful diagnostic data like this is worth it.  Hm...
```Python
In [1]: IonizationStates({'H': [0.1, 0.9], 'He': [0.9, 0.099, 0.001]})
Out[1]: 
IonizationStates instance for: H, He

  H  0+: 0.100
  H  1+: 0.900

  He  0+: 0.900
  He  1+: 0.099
```
namurphy: > Re: repr, could that perhaps be put into another custom method?

Yeah, that's what I was starting to think too, so let's go with this approach.  Thanks!
namurphy: The docstrings still need to describe the consequences of dunder methods and it'd be good to fill in a few more gaps in text coverage, but I am ready for this to be merged.  After then, I can create separate issues for the remaining tasks.
StanczakDominik: I'd just do one final round of review and likely merge this today - but what of the commit messages? I could try to trim some of it down to something more manageable. Is that okay?
namurphy: I was thinking that it'd be easiest if I did the squash merge after you approve it or I make any final changes.  That way I can remove commits that are now extraneous and maybe add some context to the remaining commits.  Does that sound good?
namurphy: Thank you, @StanczakDominik!  I am very ready for this to be merged. :joy_cat:  I am working on the summary commit message now.  

A final thought is that some of the tests could likely be refactored into fixtures.  I did not do this because I do not yet know how to write fixtures, but that is fine to leave for the future.


StanczakDominik: That sounds like a possible future low priority issue. :laughing: 

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/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",


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