Repository: PlasmaPy/PlasmaPy
Pull request #851 (2020-06-10 23:44:57)

Ionization
<!--
Thank you for contributing to PlasmaPy! Here's a bunch of pointers to
make things easier for all of us:

* If this PR will solve an issue tracked by GitHub, then please add
  "Closes #42" so the issue automatically closes once this pull request
  is merged.  In this example, issue #42 would have been closed.
  If your PR will not completely solve the issue, then please still reference
  still reference the issue like "issue #42".  (For more info see [GitHub's primer
  on linking issues](https://help.github.com/en/github/managing-your-work-on-github/linking-a-pull-request-to-an-issue#linking-a-pull-request-to-an-issue-using-a-keyword).)

* Remember to add some description of your changes in this text box.

* If your pull request is not yet ready for review - either because
  you're just looking for feedback on a change or because you're not
  perfectly satisfied with your change - submit it as a draft pull request.
  Remember to change its' status once it's ready.

* If this is your first contribution, please add your name to the author
  list in `docs/about/credits.rst`.

* Feel free to chat with other developers on our Matrix channel at:
   https://riot.im/app/#/room/#plasmapy:openastronomy.org

* We have a developer's guide to help answer some of your questions.
  http://docs.plasmapy.org/en/latest/development/index.html

Many thanks in advance for following these pointers and for being willing to contribute!

When submitting a pull request, please ensure that you can (eventually,
sometime before it is merged) check the following basic requirements:

-->

- [x] I have added a changelog entry for this pull request.

<!--

In short: A changelog entry is a short description of your PR's changes.
Each entry is written in a `<PULL REQUEST>.<TYPE>.rst` file and stored in
the `changelog` directory,  where `<PULL REQUEST>` is a pull request
number and `<TYPE>` is one of:

* `breaking`: A change which requires users to change code and is not backwards compatible. (Not to be used for removal of deprecated features.)
* `feature`: New user facing features and any new behavior.
* `bugfix`: Fixes a reported bug.
* `doc`: Documentation addition or improvement, like rewording an entire session or adding missing docs.
* `removal`: Feature deprecation and/or feature removal.
* `trivial`: A change which has no user facing effect or is tiny change.

A PR number is generated after you successfully submit a new PR, so the changelog
file can only be added after you open the PR.```

For more information, see:
https://github.com/PlasmaPy/PlasmaPy/blob/master/changelog/README.rst

--->

- [x] If adding new functionality, I have added tests and
      docstrings.

<!--
(Tests pop up at the bottom, in the checks box.)
-->

- [x] I have fixed any newly failing tests.

<!--
(If you're unsure why they're failing, ask!)
-->
This adds the `ionization` module to the formulary as well as Z_bal from the Saha equation

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

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

```diff
@@            Coverage Diff             @@
##           master     #851      +/-   ##
==========================================
- Coverage   95.94%   95.82%   -0.12%     
==========================================
  Files          56       57       +1     
  Lines        5126     5100      -26     
==========================================
- Hits         4918     4887      -31     
- Misses        208      213       +5     
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/851?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/formulary/dimensionless.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/851/diff?src=pr&el=tree#diff-cGxhc21hcHkvZm9ybXVsYXJ5L2RpbWVuc2lvbmxlc3MucHk=) | `100.00% <ø> (ø)` | |
| [plasmapy/formulary/parameters.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/851/diff?src=pr&el=tree#diff-cGxhc21hcHkvZm9ybXVsYXJ5L3BhcmFtZXRlcnMucHk=) | `98.88% <ø> (ø)` | |
| [plasmapy/formulary/ionization.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/851/diff?src=pr&el=tree#diff-cGxhc21hcHkvZm9ybXVsYXJ5L2lvbml6YXRpb24ucHk=) | `54.54% <54.54%> (ø)` | |
| [plasmapy/classes/plasma\_base.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/851/diff?src=pr&el=tree#diff-cGxhc21hcHkvY2xhc3Nlcy9wbGFzbWFfYmFzZS5weQ==) | `67.74% <0.00%> (-4.49%)` | :arrow_down: |
| [plasmapy/simulation/abstractions.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/851/diff?src=pr&el=tree#diff-cGxhc21hcHkvc2ltdWxhdGlvbi9hYnN0cmFjdGlvbnMucHk=) | `75.00% <0.00%> (-2.78%)` | :arrow_down: |
| [plasmapy/classes/sources/openpmd\_hdf5.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/851/diff?src=pr&el=tree#diff-cGxhc21hcHkvY2xhc3Nlcy9zb3VyY2VzL29wZW5wbWRfaGRmNS5weQ==) | `78.40% <0.00%> (-1.17%)` | :arrow_down: |
| [plasmapy/simulation/particletracker.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/851/diff?src=pr&el=tree#diff-cGxhc21hcHkvc2ltdWxhdGlvbi9wYXJ0aWNsZXRyYWNrZXIucHk=) | `98.27% <0.00%> (-0.03%)` | :arrow_down: |
| [plasmapy/particles/ionization\_state.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/851/diff?src=pr&el=tree#diff-cGxhc21hcHkvcGFydGljbGVzL2lvbml6YXRpb25fc3RhdGUucHk=) | `92.93% <0.00%> (-0.03%)` | :arrow_down: |
| [plasmapy/utils/pytest\_helpers/pytest\_helpers.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/851/diff?src=pr&el=tree#diff-cGxhc21hcHkvdXRpbHMvcHl0ZXN0X2hlbHBlcnMvcHl0ZXN0X2hlbHBlcnMucHk=) | `93.38% <0.00%> (-0.03%)` | :arrow_down: |
| [plasmapy/particles/particle\_class.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/851/diff?src=pr&el=tree#diff-cGxhc21hcHkvcGFydGljbGVzL3BhcnRpY2xlX2NsYXNzLnB5) | `95.73% <0.00%> (-0.01%)` | :arrow_down: |
| ... and [8 more](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/851/diff?src=pr&el=tree-more) | |

------

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

StanczakDominik: Remember to commit/push the alias change from today, and this will probably be a quick merge! :)

Patch:
diff --git a/changelog/831.feature.rst b/changelog/831.feature.rst
new file mode 100644
index 00000000..6ae1b3b9
--- /dev/null
+++ b/changelog/831.feature.rst
@@ -0,0 +1,5 @@
+Added a new diagnostics module `thomson` containing a function
+`spectral_density` that calculates Thomson scattering spectra for
+Maxwellian plasmas in both the collective and non-collective regimes. As
+a followup to PR #835, set the minimal required Numpy version to 1.18.1 to
+finally fix unit dropping bugs.
diff --git a/changelog/850.feature.rst b/changelog/850.feature.rst
new file mode 100644
index 00000000..23bd21f6
--- /dev/null
+++ b/changelog/850.feature.rst
@@ -0,0 +1,1 @@
+Revised parameters.thermal_speed to support 1D and 2D distributions as well as 3D, and added an example notebook for this function.
diff --git a/changelog/851.feature.rst b/changelog/851.feature.rst
new file mode 100644
index 00000000..219fe34d
--- /dev/null
+++ b/changelog/851.feature.rst
@@ -0,0 +1,2 @@
+Create `plasmapy/formulary/ionization.py`
+Create :func:`~plasmapy.formulary.ionization.Z_bal` function.
diff --git a/docs/about/credits.rst b/docs/about/credits.rst
index 886b091c..604b47f6 100644
--- a/docs/about/credits.rst
+++ b/docs/about/credits.rst
@@ -42,6 +42,7 @@ in parentheses are `ORCID author identifiers <https://orcid.org>`__.
 * `Graham Goudeau <https://github.com/GrahamGoudeau>`__
 * `Silvina Guidoni <https://www.american.edu/cas/faculty/guidoni.cfm>`__
 * `Colby Haggerty <https://github.com/colbych>`__
+* `Peter Heuer <https://github.com/pheuer>`__ (`0000-0001-5050-6606 <https://orcid.org/0000-0001-5050-6606>`__)
 * `Julien Hillairet <https://github.com/jhillairet>`__ (`0000-0002-1073-6383 <https://orcid.org/0000-0002-1073-6383>`__)
 * `Poh Zi How <https://github.com/pohzipohzi>`__
 * `Yi-Min Huang <https://github.com/yopology>`__ (`0000-0002-4237-2211 <https://orcid.org/0000-0002-4237-2211>`__)
diff --git a/docs/conf.py b/docs/conf.py
index 34a7fb01..e251d404 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -24,9 +24,10 @@
 from pkg_resources import parse_version
 from sphinx.application import Sphinx
 
+sys.path.insert(0, os.path.abspath(".."))
+
 from plasmapy import __version__ as release
 
-sys.path.insert(0, os.path.abspath(".."))
 
 # -- General configuration ------------------------------------------------
 
diff --git a/docs/diagnostics/index.rst b/docs/diagnostics/index.rst
index 09bbd749..d5d32128 100644
--- a/docs/diagnostics/index.rst
+++ b/docs/diagnostics/index.rst
@@ -4,14 +4,12 @@
 Plasma diagnostics (`plasmapy.diagnostics`)
 *******************************************
 
-.. currentmodule:: plasmapy.diagnostics.langmuir
+.. automodapi:: plasmapy.diagnostics
+    :no-heading:
 
-.. automodapi:: plasmapy.diagnostics.langmuir
-   :no-heading:
+.. toctree::
+   :maxdepth: 1
+   :caption: Sub-Packages and Modules
 
-|
-
-.. nbgallery::
-    :caption: Examples
-
-    /notebooks/langmuir_analysis
+   langmuir
+   thomson
diff --git a/docs/diagnostics/langmuir.rst b/docs/diagnostics/langmuir.rst
new file mode 100644
index 00000000..3d71ec16
--- /dev/null
+++ b/docs/diagnostics/langmuir.rst
@@ -0,0 +1,17 @@
+.. _langmuir:
+
+*****************
+Langmuir analysis
+*****************
+
+.. currentmodule:: plasmapy.diagnostics.langmuir
+
+.. automodapi:: plasmapy.diagnostics.langmuir
+   :no-heading:
+
+|
+
+.. nbgallery::
+    :caption: Examples
+
+    /notebooks/langmuir_analysis
diff --git a/docs/diagnostics/thomson.rst b/docs/diagnostics/thomson.rst
new file mode 100644
index 00000000..96493b3e
--- /dev/null
+++ b/docs/diagnostics/thomson.rst
@@ -0,0 +1,15 @@
+.. _thomson:
+
+******************
+Thomson scattering
+******************
+
+.. currentmodule:: plasmapy.diagnostics.thomson
+
+.. automodapi:: plasmapy.diagnostics.thomson
+   :no-heading:
+
+.. nbgallery::
+    :caption: Examples
+
+    /notebooks/thomson
diff --git a/docs/formulary/distribution.rst b/docs/formulary/distribution.rst
index c805261e..2cc487b7 100644
--- a/docs/formulary/distribution.rst
+++ b/docs/formulary/distribution.rst
@@ -13,3 +13,4 @@ Distribution functions (`plasmapy.formulary.distribution`)
     :caption: Examples
 
     /notebooks/distribution
+    /notebooks/thermal_speed
diff --git a/docs/formulary/parameters.rst b/docs/formulary/parameters.rst
index 76bd8de4..65ba803c 100644
--- a/docs/formulary/parameters.rst
+++ b/docs/formulary/parameters.rst
@@ -11,3 +11,4 @@ Plasma parameters (`plasmapy.formulary.parameters`)
     :caption: Examples
 
     /notebooks/distribution
+    /notebooks/thermal_speed
diff --git a/docs/install.rst b/docs/install.rst
index 72b776a9..4aad2653 100644
--- a/docs/install.rst
+++ b/docs/install.rst
@@ -18,7 +18,7 @@ Requirements
 PlasmaPy requires Python version 3.6 or newer.
 PlasmaPy requires the following packages for installation:
 
-- `NumPy <http://www.numpy.org/>`_ 1.16.6 or newer
+- `NumPy <http://www.numpy.org/>`_ 1.18.1 or newer
 - `SciPy <https://www.scipy.org/>`_ 1.2 or newer
 - `Astropy <http://www.astropy.org/>`_ 4.0 or newer
 - `colorama <https://pypi.org/project/colorama/>`_ 0.3 or newer
diff --git a/docs/notebooks/thermal_speed.ipynb b/docs/notebooks/thermal_speed.ipynb
new file mode 100644
index 00000000..4dc999d9
--- /dev/null
+++ b/docs/notebooks/thermal_speed.ipynb
@@ -0,0 +1,156 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Thermal Speed"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "\n",
+    "import numpy as np\n",
+    "from astropy import units as u\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "from plasmapy.formulary import  Maxwellian_speed_1D, Maxwellian_speed_2D, Maxwellian_speed_3D\n",
+    "from plasmapy.formulary.parameters import thermal_speed"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The thermal_speed function can be used to calculate the thermal velocity for a Maxwellian velocity distribution. There are three common definitions of the thermal velocity, which can be selected using the \"method\" keyword, which are defined for a 3D velocity distribution as\n",
+    "\n",
+    "\n",
+    "- 'most_probable' <br>\n",
+    "$v_{th} = \\sqrt{\\frac{2 k_B T}{m}}$\n",
+    "\n",
+    "- 'rms' <br>\n",
+    "$v_{th} = \\sqrt{\\frac{3 k_B T}{m}}$\n",
+    "\n",
+    "- 'mean_magnitude' <br>\n",
+    "$v_{th} = \\sqrt{\\frac{8 k_B T}{m\\pi}}$\n",
+    "\n",
+    "The differences between these velocities can be seen by plotitng them on a 3D Maxwellian speed distribution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "nbsphinx-thumbnail": {
+    "tooltip": "Thermal Speeds"
+    }
+    },
+   "outputs": [],
+   "source": [
+    "T = 1e5 * u.K\n",
+    "speeds = np.linspace(0, 8e6, num=600) * u.m/u.s\n",
+    "\n",
+    "pdf_3D = Maxwellian_speed_3D(speeds, T=T, particle='e-')\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(4,3))\n",
+    "\n",
+    "v_most_prob = thermal_speed(T=T, particle='e-', method='most_probable', ndim=3)\n",
+    "v_rms = thermal_speed(T=T, particle='e-', method='rms', ndim=3)\n",
+    "v_mean_magnitude = thermal_speed(T=T, particle='e-', method='mean_magnitude', ndim=3)\n",
+    "\n",
+    "ax.plot(speeds/v_rms, pdf_3D, color='black', label='Maxwellian')\n",
+    "    \n",
+    "ax.axvline(x=v_most_prob/v_rms, color='blue', label='Most Probable')\n",
+    "ax.axvline(x=v_rms/v_rms, color='green', label='RMS')\n",
+    "ax.axvline(x=v_mean_magnitude/v_rms, color='red', label='Mean Magnitude')\n",
+    "\n",
+    "ax.set_xlim(-.1, 3)\n",
+    "ax.set_ylim(0, None)\n",
+    "ax.set_title('3D')\n",
+    "ax.set_xlabel(\"|v|/|v$_{rms}|$\")\n",
+    "ax.set_ylabel(\"f(|v|)\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Similar speeds are defined for 1D and 2D distributions. The differences between these definitions can be illustrated by plotting them on their respective Maxwellian speed distributions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "pdf_1D = Maxwellian_speed_1D(speeds, T=T, particle='e-')\n",
+    "pdf_2D = Maxwellian_speed_2D(speeds, T=T, particle='e-')\n",
+    "\n",
+    "dim = [1,2,3]\n",
+    "pdfs = [pdf_1D, pdf_2D, pdf_3D]\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "fig, ax = plt.subplots(ncols=3, figsize=(10,3))\n",
+    "\n",
+    "for n, pdf in enumerate(pdfs):\n",
+    "    ndim = n+1\n",
+    "    v_most_prob = thermal_speed(T=T, particle='e-', method='most_probable', ndim=ndim)\n",
+    "    v_rms = thermal_speed(T=T, particle='e-', method='rms', ndim=ndim)\n",
+    "    v_mean_magnitude = thermal_speed(T=T, particle='e-', method='mean_magnitude', ndim=ndim)\n",
+    "    \n",
+    "    ax[n].plot(speeds/v_rms, pdf, color='black', label='Maxwellian')\n",
+    "    \n",
+    "    ax[n].axvline(x=v_most_prob/v_rms, color='blue', label='Most Probable')\n",
+    "    ax[n].axvline(x=v_rms/v_rms, color='green', label='RMS')\n",
+    "    ax[n].axvline(x=v_mean_magnitude/v_rms, color='red', label='Mean Magnitude')\n",
+    "    \n",
+    "    ax[n].set_xlim(-.1, 3)\n",
+    "    ax[n].set_ylim(0, None)\n",
+    "    ax[n].set_title('{:d}D'.format(ndim))\n",
+    "    ax[n].set_xlabel(\"|v|/|v$_{rms}|$\")\n",
+    "    ax[n].set_ylabel(\"f(|v|)\")\n",
+    "\n",
+    "\n",
+    "ax[2].legend(bbox_to_anchor=(1.9, .8), loc='upper right')\n",
+    "\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/docs/notebooks/thomson.ipynb b/docs/notebooks/thomson.ipynb
new file mode 100644
index 00000000..a26720cb
--- /dev/null
+++ b/docs/notebooks/thomson.ipynb
@@ -0,0 +1,348 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Thomson Scattering: Spectral Density"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "[thomson]: ../diagnostics/thomson.rst\n",
+    "[spectral-density]: ../api/plasmapy.diagnostics.thomson.spectral_density.rst#spectral-density\n",
+    "[sheffield]: https://www.sciencedirect.com/book/9780123748775/plasma-scattering-of-electromagnetic-radiation\n",
+    "\n",
+    "The [thomson.spectral_density][spectral-density] function calculates the [spectral density function S(k,w)][sheffield], which is one of several terms that determine the scattered power spectrum for the Thomson scattering of a probe laser beam by a plasma. In particular, this function calculates $S(k,w)$ for the case of a plasma consisting of one or more ion species and a neutralizing electron fluid under the assumption that all of the ion species and the electron fluid have Maxwellian velocity distribution functions. In this regime, the spectral density is given by the equation:\n",
+    "\n",
+    "\\begin{equation}\n",
+    "S(k,\\omega) = \\frac{2\\pi}{k} \\bigg |1 - \\frac{\\chi_e}{\\epsilon} \\bigg |^2 f_{e0}\\bigg ( \\frac{\\omega}{k} \\bigg ) + \\sum_i \\frac{2\\pi Z_i}{k} \\bigg | \\frac{\\chi_e}{\\epsilon} \\bigg |^2 f_{i0, i} \\bigg ( \\frac{\\omega}{k} \\bigg )\n",
+    "\\end{equation}\n",
+    "\n",
+    "where $\\chi_e$ is the electron component susceptibility of the plasma and $\\epsilon = 1 + \\chi_e + \\sum_i \\chi_i$ is the total plasma dielectric function (with $\\chi_i$ being the ion component of the susceptibility), $Z_i$ is the charge of each ion, $k$ is the scattering wavenumber, $\\omega$ is the scattering frequency, and the functions $f_{e0}$ and $f_{i0,i}$ are the Maxwellian velocity distributions for the electrons and ion species respectively.\n",
+    "\n",
+    "Thomson scattering can be either non-collective (the scattered spectrum is a linear sum of the light scattered by individual particles) or collective (the scattered spectrum is dominated by scattering off of collective plasma waves). The [thomson.spectral_density][spectral-density] function can be used in both cases. These regimes are delineated by the dimensionless constant $\\alpha$:\n",
+    "\n",
+    "\\begin{equation}\n",
+    "\\alpha = \\frac{1}{k \\lambda_{De}}\n",
+    "\\end{equation}\n",
+    "\n",
+    "where $\\lambda_{De}$ is the Debye length. $\\alpha > 1$ corresponds to collective scattering, while $\\alpha < 1$ corresponds to non-collective scattering. Depending on which of these regimes applies, fitting the scattered spectrum can provide the electron (and sometimes ion) density and temperature. Doppler shifting of the spectrum can also provide a measurement of the drift velocity of each plasma species.\n",
+    "\n",
+    "For a detailed explanation of the underlying physics (and derivations of these expressions), see [\"Plasma Scattering of Electromagnetic Radiation\" by Sheffield et al.][sheffield]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "\n",
+    "import astropy.units as u\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import warnings\n",
+    "\n",
+    "from plasmapy.diagnostics import thomson\n",
+    "from plasmapy.utils.exceptions import ImplicitUnitConversionWarning"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Construct parameters that define the Thomson diagnostic setup, the probing beam and scattering collection.  These parameters will be used for all examples."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# The probe wavelength can in theory be anything, but in practice integer frequency multiples of the Nd:YAG wavelength\n",
+    "# 1064 nm are used (532 corresponds to a frequency-doubled probe beam from such a laser).\n",
+    "probe_wavelength = 532*u.nm\n",
+    "\n",
+    "# Array of wavelengths over which to calcualte the spectral distribution\n",
+    "wavelengths = np.arange(probe_wavelength.value-60, probe_wavelength.value+60, 0.01)*u.nm\n",
+    "\n",
+    "# The scattering geometry is defined by unit vectors for the orientation of the probe laser beam (probe_n) and\n",
+    "# the path from the scattering volume (where the measurement is made) to the detector (scatter_n).\n",
+    "# These can be setup for any experimental geometry.\n",
+    "probe_vec = np.array([1, 0, 0])\n",
+    "scattering_angle = np.deg2rad(63)\n",
+    "scatter_vec = np.array([np.cos(scattering_angle), np.sin(scattering_angle), 0])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In order to calcluate the scattered spectrum, we must also include some information about the plasma. For this plot we'll allow the ``fract``, ``ion_species``, ``fluid_vel``, and ``ion_vel`` keywords to keep their default values, describing a single-species H+ plasma at rest in the laboratory frame. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "nbsphinx-thumbnail": {
+     "tooltip": "Thomson Scattering"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "ne = 2e17*u.cm**-3\n",
+    "Te = 12*u.eV\n",
+    "Ti = 10*u.eV\n",
+    "\n",
+    "# This warning filter catches an ImplicitUnitConversionWarning that results from specifying\n",
+    "# temperatures in eV instead of Kelvin. \n",
+    "with warnings.catch_warnings():\n",
+    "    warnings.simplefilter(\"ignore\", ImplicitUnitConversionWarning)\n",
+    "    \n",
+    "    alpha, Skw = thomson.spectral_density(wavelengths, probe_wavelength,\n",
+    "                                          ne, Te, Ti, probe_vec=probe_vec, \n",
+    "                                          scatter_vec=scatter_vec)\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(wavelengths, Skw, lw=2)\n",
+    "ax.set_xlim(probe_wavelength.value-10, probe_wavelength.value+10)\n",
+    "ax.set_ylim(0, 1e-13)\n",
+    "ax.set_xlabel('$\\lambda$ (nm)')\n",
+    "ax.set_ylabel('S(k,w)')\n",
+    "ax.set_title('Thomson Scattering Spectral Density')\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Example Cases in Different Scattering Regimes\n",
+    "\n",
+    "We will now consider several example cases in different scattering regimes. In order to facilitate this, we'll set up each example as a dictionary of plasma parameters:"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A single-species, stationary hydrogen plasma with a density and temperature that results in a scattering spectrum dominated by scattering off of single electrons."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "non_collective = {\n",
+    "    'name': 'Non-Collective Regime', \n",
+    "    'ne': 5e15*u.cm**-3, \n",
+    "    'Te': 40*u.eV, \n",
+    "    'Ti': np.array([10])*u.eV,\n",
+    "    'fract': np.array([1]), \n",
+    "    'ion_species': ['H+'],\n",
+    "    'fluid_vel': np.array([0, 0, 0])*u.km/u.s, \n",
+    "    'ion_vel': np.array([[0, 0, 0]])*u.km/u.s,\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A single-species, stationary hydrogen plasma with a density and temperature that result in weakly collective scattering (scattering paramter $\\alpha$ approaching 1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "weakly_collective = {\n",
+    "    'name': 'Weakly Collective Regime',\n",
+    "    'ne': 2e17*u.cm**-3,\n",
+    "    'Te': 20*u.eV,\n",
+    "    'Ti': np.array([10])*u.eV,\n",
+    "    'fract': np.array([1]),\n",
+    "    'ion_species': ['H+'],\n",
+    "    'fluid_vel': np.array([0, 0, 0])*u.km/u.s,\n",
+    "    'ion_vel': np.array([[0, 0, 0]])*u.km/u.s,\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A single-species, stationary hydrogen plasma with a density and temperature that result in a spectrum dominated by multi-particle scattering, including scattering off of ions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "collective = {\n",
+    "    'name': 'Collective Regime',\n",
+    "    'ne': 5e17*u.cm**-3,\n",
+    "    'Te': 10*u.eV,\n",
+    "    'Ti': np.array([4])*u.eV,\n",
+    "    'fract': np.array([1]),\n",
+    "    'ion_species': ['H+'],\n",
+    "    'fluid_vel': np.array([0, 0, 0])*u.km/u.s,\n",
+    "    'ion_vel': np.array([[0, 0, 0]])*u.km/u.s,\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A case identical to the collective example above, except that now the electron fluid has a substantial drift velocity parallel to the probe laser and the ions have a drift (relative to the electrons) at an angle."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "drifts = {\n",
+    "    'name': 'Drift Velocities',\n",
+    "    'ne': 5e17*u.cm**-3,\n",
+    "    'Te': 10*u.eV,\n",
+    "    'Ti': np.array([10])*u.eV,\n",
+    "    'fract': np.array([1]),\n",
+    "    'ion_species': ['H+'],\n",
+    "    'fluid_vel': np.array([700, 0, 0])*u.km/u.s,\n",
+    "    'ion_vel': np.array([[-600, -100, 0]])*u.km/u.s,\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A case identical to the collective example, except that now the plasma consists 25% He+1 and 75% C+5"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "two_species = {\n",
+    "    'name': 'Two Ion Species',\n",
+    "    'ne': 5e17*u.cm**-3,\n",
+    "    'Te': 10*u.eV,\n",
+    "    'Ti': np.array([10, 50])*u.eV,\n",
+    "    'fract': np.array([.25, .75]),\n",
+    "    'ion_species': ['He-4 1+', 'C-12 5+'],\n",
+    "    'fluid_vel': np.array([0, 0, 0])*u.km/u.s,\n",
+    "    'ion_vel': np.array([[0, 0, 0],[0, 0, 0]])*u.km/u.s,\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "examples = [non_collective, weakly_collective, collective, drifts, two_species]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For each example, plot the the spectral distribution function over a large range to show the broad electron scattering feature (top row) and a narrow range around the probe wavelength to show the ion scattering feature (bottom row)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(ncols=len(examples), nrows=2, figsize=[25,10])\n",
+    "fig.subplots_adjust( wspace=0.4, hspace=0.4)\n",
+    "\n",
+    "lbls = 'abcdefg'\n",
+    "\n",
+    "for i, x in enumerate(examples):\n",
+    "    \n",
+    "    with warnings.catch_warnings():\n",
+    "        warnings.simplefilter(\"ignore\", ImplicitUnitConversionWarning)\n",
+    "        \n",
+    "        alpha, Skw = thomson.spectral_density(wavelengths, probe_wavelength,\n",
+    "                         x['ne'], x['Te'], x['Ti'], fract=x['fract'],\n",
+    "                         ion_species=x['ion_species'], fluid_vel=x['fluid_vel'],\n",
+    "                         probe_vec=probe_vec, scatter_vec=scatter_vec)\n",
+    "    \n",
+    "    ax[0][i].axvline(x=probe_wavelength.value, color='red')  # Mark the probe wavelength\n",
+    "    ax[0][i].plot(wavelengths, Skw)\n",
+    "    ax[0][i].set_xlim(probe_wavelength.value-15, probe_wavelength.value+15)\n",
+    "    ax[0][i].set_ylim(0, 1e-13)\n",
+    "    ax[0][i].set_xlabel('$\\lambda$ (nm)')\n",
+    "    \n",
+    "    ax[0][i].set_title( lbls[i] + ') ' + x['name'] + '\\n$\\\\alpha$={:.4f}'.format(alpha))\n",
+    "\n",
+    "    \n",
+    "    ax[1][i].axvline(x=probe_wavelength.value, color='red')  # Mark the probe wavelength\n",
+    "    ax[1][i].plot(wavelengths, Skw)\n",
+    "    ax[1][i].set_xlim(probe_wavelength.value-1, probe_wavelength.value+1)\n",
+    "    ax[1][i].set_ylim(0, 1.1*np.max(Skw.value))\n",
+    "    ax[1][i].set_xlabel('$\\lambda$ (nm)')\n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plots of the spectral density function (Skw) which determines the amount of light scattered into different wavelengths.\n",
+    "\n",
+    "a. In the non-collective regime only the electron feature is visible.\n",
+    "\n",
+    "b. In the weakly collective regime (alpha approaches 1) an ion feature starts to appear and the electron feature is distorted\n",
+    "\n",
+    "c. In the collective regime both features split into two peaks, corresponding to scattering off of forward and backwards propagating plasma oscillations. \n",
+    "\n",
+    "d. The introduction of drift velocities introduces several Doppler shifts in the calculations, resulting in a shifted spectrum.\n",
+    "\n",
+    "e. Changing the ion species present changes the resulting ion feature.\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/plasmapy/diagnostics/__init__.py b/plasmapy/diagnostics/__init__.py
index 9d8b717c..79305478 100644
--- a/plasmapy/diagnostics/__init__.py
+++ b/plasmapy/diagnostics/__init__.py
@@ -1,3 +1,7 @@
-"""The diagnostics subpackage contains tools for experimental research.
-Currently, we have functionality for analyzing data from Langmuir probes.
 """
+The diagnostics subpackage contains functionality for defining diagnostic parameters
+and processing data collected by diagnostics, synthetic or experimental.
+"""
+__all__ = ["langmuir", "thomson"]
+
+from plasmapy.diagnostics import langmuir, thomson
diff --git a/plasmapy/diagnostics/thomson.py b/plasmapy/diagnostics/thomson.py
new file mode 100644
index 00000000..41d29e6c
--- /dev/null
+++ b/plasmapy/diagnostics/thomson.py
@@ -0,0 +1,277 @@
+"""
+Defines the Thomson scattering analysis module as
+part of the diagnostics package.
+"""
+
+__all__ = [
+    "spectral_density",
+]
+
+import astropy.constants as const
+import astropy.units as u
+import numpy as np
+
+from plasmapy.formulary.parameters import plasma_frequency, thermal_speed
+from plasmapy.formulary.dielectric import permittivity_1D_Maxwellian
+from plasmapy.particles import Particle
+from plasmapy.utils.decorators import validate_quantities
+from typing import List, Tuple, Union
+
+
+# TODO: interface for inputting a multi-species configuration could be
+# simplified using the plasmapy.classes.plasma_base class if that class
+# included ion and electron drift velocities and information about the ion
+# atomic species.
+
+
+@validate_quantities(
+    wavelengths={"can_be_negative": False},
+    probe_wavelength={"can_be_negative": False},
+    ne={"can_be_negative": False},
+    Te={"can_be_negative": False, "equivalencies": u.temperature_energy()},
+    Ti={"can_be_negative": False, "equivalencies": u.temperature_energy()},
+)
+def spectral_density(
+    wavelengths: u.nm,
+    probe_wavelength: u.nm,
+    ne: u.m ** -3,
+    Te: u.K,
+    Ti: u.K,
+    fract: np.ndarray = np.ones(1),
+    ion_species: Union[str, List[str], Particle, List[Particle]] = "H+",
+    fluid_vel: u.m / u.s = np.zeros(3) * u.m / u.s,
+    ion_vel: u.m / u.s = None,
+    probe_vec=np.array([1, 0, 0]),
+    scatter_vec=np.array([0, 1, 0]),
+) -> Tuple[Union[np.floating, np.ndarray], np.ndarray]:
+    r"""
+    Calculate the spectral density function for Thomson scattering of a
+    probe laser beam by a multi-species Maxwellian plasma.
+
+    This function calculates the spectral density function for Thomson
+    scattering of a probe laser beam by a plasma consisting of one or more ion
+    species and a neutralizing electron fluid:
+
+    .. math::
+        S(k,\omega) = 1 - \frac{2\pi}{k}
+        \bigg |1 - \frac{\chi_e}{\epsilon} \bigg |^2
+        f_{e0} \bigg (\frac{\omega}{k} \bigg ) +
+        \sum_i \frac{2\pi Z_i}{k}
+        \bigg |\frac{\chi_e}{\epsilon} \bigg |^2 f_{i0,i}
+        \bigg ( \frac{\omega}{k} \bigg )
+
+    where :math:`\chi_e` is the electron component susceptibility of the
+    plasma and :math:`\epsilon = 1 + \chi_e + \sum_i \chi_i` is the total
+    plasma dielectric  function (with :math:`\chi_i` being the ion component
+    of the susceptibility), :math:`Z_i` is the charge of each ion, :math:`k`
+    is the scattering wavenumber, :math:`\omega` is the scattering frequency,
+    and :math:`f_{e0}` and :math:`f_{i0,i}` are the electron and ion velocity
+    distribution functions respectively. In this function the electron and ion
+    velocity distribution functions are assumed to be Maxwellian, making this
+    function equivalent to Eq. 3.4.6 in `Sheffield`_.
+
+    Parameters
+    ----------
+
+    wavelengths : `~astropy.units.Quantity`
+        Array of wavelengths over which the spectral density function
+        will be calculated. (convertible to nm)
+
+    probe_wavelength : `~astropy.units.Quantity`
+        Wavelength of the probe laser. (convertible to nm)
+
+    ne : `~astropy.units.Quantity`
+        Mean (0th order) electron density of all plasma components combined.
+        (convertible to cm^-3.)
+
+    Te : `~astropy.units.Quantity`
+        Temperature of the electron component. (in K or convertible to eV)
+
+    Ti : `~astropy.units.Quantity`, shape (N, )
+        Temperature of each ion component. Shape (N, ) must be equal to the
+        number of ion components N. (in K or convertible to eV)
+
+    fract : array_like, shape (N, ), optional
+        An array-like object where each element represents the fraction (or ratio)
+        of the ion component number density to the total ion number density.
+        Must sum to 1.0. Default is a single ion species.
+
+    ion_species : str or `~plasmapy.particles.Particle`, shape (N, ), optional
+        A list or single instance of `~plasmapy.particles.Particle`, or strings
+        convertible to `~plasmapy.particles.Particle`. Default is `'H+'`
+        corresponding to a single species of hydrogen ions.
+
+    fluid_vel : `~astropy.units.Quantity`, shape (3, ), optional
+        Electron fluid velocity in the rest frame. (convertible to m/s)
+        Defaults to a stationary plasma [0, 0, 0] m/s.
+
+    ion_vel : `~astropy.units.Quantity`, shape (N, 3), optional
+        Velocity vectors for each ion population relative to the
+        fluid velocity. (convertible to m/s) Defaults zero drift
+        for all specified ion species.
+
+    probe_vec : float `~numpy.ndarray`, shape (3, )
+        Unit vector in the direction of the probe laser. Defaults to
+        [1, 0, 0].
+
+    scatter_vec : float `~numpy.ndarray`, shape (3, )
+        Unit vector pointing from the scattering volume to the detector.
+        Defaults to [0, 1, 0] which, along with the default `probe_vec`,
+        corresponds to a 90 degree scattering angle geometry.
+
+    Returns
+    -------
+    alpha : float
+        Mean scattering parameter, where `alpha` > 1 corresponds to collective
+        scattering and `alpha` < 1 indicates non-collective scattering.
+
+    Skw : `~astropy.units.Quantity`
+        Computed spectral density function over the input `wavelengths` array
+        with units of s/rad.
+
+    Notes
+    -----
+
+    For details, see "Plasma Scattering of Electromagnetic Radiation" by
+    Sheffield et al. `ISBN 978\\-0123748775`_. This code is a modified version
+    of the program described therein.
+
+    For a concise summary of the relevant physics, see Chapter 5 of Derek
+    Schaeffer's thesis, DOI: `10.5281/zenodo.3766933`_.
+
+    .. _`ISBN 978\\-0123748775`: https://www.sciencedirect.com/book/9780123748775/plasma-scattering-of-electromagnetic-radiation
+    .. _`10.5281/zenodo.3766933`: https://doi.org/10.5281/zenodo.3766933
+    .. _`Sheffield`: https://doi.org/10.1016/B978-0-12-374877-5.00003-8
+    """
+    fract = np.asarray(fract, dtype=np.float64)
+
+    # If ion drift velocity is not specified, create an array corresponding
+    # to zero drift
+    if ion_vel is None:
+        ion_vel = np.zeros([fract.size, 3]) * u.m / u.s
+
+    # Condition ion_species
+    if isinstance(ion_species, (str, Particle)):
+        ion_species = [ion_species]
+    if len(ion_species) == 0:
+        raise ValueError("At least one ion species needs to be defined.")
+    for ii, ion in enumerate(ion_species):
+        if isinstance(ion, Particle):
+            continue
+        ion_species[ii] = Particle(ion)
+
+    # Condition Ti
+    if Ti.size == 1:
+        # If a single quantity is given, put it in an array so it's iterable
+        # If Ti.size != len(ion_species), assume same temp. for all species
+        Ti = [Ti.value] * len(ion_species) * Ti.unit
+    elif Ti.size != len(ion_species):
+        raise ValueError(
+            f"Got {Ti.size} ion temperatures and expected " f"{len(ion_species)}."
+        )
+
+    if (
+        (len(ion_species) != fract.size)
+        or (ion_vel.shape[0] != fract.size)
+        or (Ti.size != fract.size)
+    ):
+        raise ValueError(
+            "Inconsistent number of species in fract, "
+            "ion_species, Ti, and/or ion_vel."
+        )
+
+    # Ensure unit vectors are normalized
+    probe_vec = probe_vec / np.linalg.norm(probe_vec)
+    scatter_vec = scatter_vec / np.linalg.norm(scatter_vec)
+
+    # Define some constants
+    C = const.c.si  # speed of light
+
+    # Calculate plasma parameters
+    vTe = thermal_speed(Te, particle="e-")
+    vTi, ion_z = [], []
+    for T, ion in zip(Ti, ion_species):
+        vTi.append(thermal_speed(T, particle=ion).value)
+        ion_z.append(ion.integer_charge * u.dimensionless_unscaled)
+    vTi = vTi * vTe.unit
+    zbar = np.sum(fract * ion_z)
+    ni = fract * ne / zbar  # ne/zbar = sum(ni)
+    wpe = plasma_frequency(n=ne, particle="e-")
+
+    # Compute the ion velocity in the rest frame
+    ion_vel = fluid_vel + ion_vel
+
+    # Convert wavelengths to angular frequencies (electromagnetic waves, so
+    # phase speed is c)
+    ws = (2 * np.pi * u.rad * C / wavelengths).to(u.rad / u.s)
+    wl = (2 * np.pi * u.rad * C / probe_wavelength).to(u.rad / u.s)
+
+    # Compute the frequency shift (required by energy conservation)
+    w = ws - wl
+
+    # Compute the wavenumbers in the plasma
+    # See Sheffield Sec. 1.8.1 and Eqs. 5.4.1 and 5.4.2
+    ks = np.sqrt(ws ** 2 - wpe ** 2) / C
+    kl = np.sqrt(wl ** 2 - wpe ** 2) / C
+
+    # Compute the wavenumber shift (required by momentum conservation)
+    scattering_angle = np.arccos(np.dot(probe_vec, scatter_vec))
+    # Eq. 1.7.10 in Sheffield
+    k = np.sqrt(ks ** 2 + kl ** 2 - 2 * ks * kl * np.cos(scattering_angle))
+    # Normal vector along k
+    k_vec = (scatter_vec - probe_vec) * u.dimensionless_unscaled
+
+    # Compute Doppler-shifted frequencies for both the ions and electrons
+    # Matmul is simultaneously conducting dot product over all wavelengths
+    # and ion components
+    w_e = w - k * np.dot(fluid_vel, k_vec)
+    w_i = w - np.matmul(ion_vel, np.outer(k, k_vec).T)
+
+    # Compute the scattering parameter alpha
+    # expressed here using the fact that v_th/w_p = root(2) * Debye length
+    alpha = np.sqrt(2) * wpe / (k * vTe)
+
+    # Calculate the normalized phase velocities (Sec. 3.4.2 in Sheffield)
+    xe = (w_e / (k * vTe)).to(u.dimensionless_unscaled)
+    xi = (np.outer(1 / vTi, 1 / k) * w_i).to(u.dimensionless_unscaled)
+
+    # Calculate the susceptibilities
+    chiE = permittivity_1D_Maxwellian(w_e, k, Te, ne, "e-")
+
+    # Treatment of multiple species is an extension of the discussion in
+    # Sheffield Sec. 5.1
+    chiI = np.zeros([fract.size, w.size], dtype=np.complex128)
+    for i, ion in enumerate(ion_species):
+        chiI[i, :] = permittivity_1D_Maxwellian(
+            w_i[i, :], k, Ti[i], ni[i], ion, z_mean=ion_z[i]
+        )
+
+    # Calculate the longitudinal dielectric function
+    epsilon = 1 + chiE + np.sum(chiI, axis=0)
+
+    # Calculate the contributions to the spectral density function
+    econtr = (
+        2
+        * np.sqrt(np.pi)
+        / k
+        / vTe
+        * np.power(np.abs(1 - chiE / epsilon), 2)
+        * np.exp(-(xe ** 2))
+    )
+
+    icontr = np.zeros([fract.size, w.size], dtype=np.complex128) * u.s / u.rad
+    for m in range(fract.size):
+        icontr[m, :] = (
+            2
+            * np.sqrt(np.pi)
+            * ion_z[m]
+            / k
+            / vTi[m]
+            * np.power(np.abs(chiE / epsilon), 2)
+            * np.exp(-xi[m, :] ** 2)
+        )
+
+    # Recast as real: imaginary part is already zero
+    Skw = np.real(econtr + np.sum(icontr, axis=0))
+
+    return np.mean(alpha), Skw
diff --git a/plasmapy/formulary/ionization.py b/plasmapy/formulary/ionization.py
new file mode 100644
index 00000000..5909b861
--- /dev/null
+++ b/plasmapy/formulary/ionization.py
@@ -0,0 +1,90 @@
+""" This module gathers functions relating to ionization states and the properties thereof.
+"""
+__all__ = ["ionization_balance"]
+
+from astropy import units as u
+from plasmapy.utils.decorators import (
+    angular_freq_to_hz,
+    check_relativistic,
+    validate_quantities,
+)
+from astropy.constants import c, a0, k_B
+from numpy import pi, exp, sqrt, log
+
+
+@validate_quantities(
+    T_e={"can_be_negative": False, "equivalencies": u.temperature_energy()}
+)
+def ionization_balance(n: u.m ** -3, T_e: u.K) -> u.dimensionless_unscaled:
+    r"""
+    Z_bal is the estimate average ionization level of a plasma in thermal equilibrium
+    that balances the number density of ions in two different ionization states.
+    Z_bal is derived from the Saha equation with the assumptions that
+    the atoms are of a single species,
+    are either hydrogenic or completely ionized,
+    and that there is a balance between ionization and recombination,
+    meaning that the number of atoms in either state are equal.
+    The Saha equation and therefore Z_bal are more accurate when the plasma
+    is at a high density and temperature.
+
+    .. math::
+
+        Z\_bal = \sqrt{\frac{k_B T_e}{E_H}} \sqrt{\ln{\frac{1}{4 n a_{0}^3}
+        (\frac{k_B T_e}{\pi E_H})^{3/2}}} - \frac{1}{2}
+
+    Where :math:`k_B` is the Boltzmann constant,
+    :math:`a_0` is the Bohr radius, and
+    :math:`E_H` is the ionization energy of Hydrogen
+
+    Parameters
+    ----------
+    T_e : `~astropy.units.Quantity`
+        The electron temperature.
+
+    n : `~astropy.units.Quantity`
+        The electron number density of the plasma.
+
+    Warns
+    -----
+    ~astropy.units.UnitsWarning
+        If units are not provided, SI units are assumed.
+
+    Raises
+    ------
+    TypeError
+        The `T_e` or `n` is not a `~astropy.units.Quantity` and cannot be
+        converted into a ~astropy.units.Quantity.
+
+    ~astropy.units.UnitConversionError
+        If the `T_e` or `n` not in appropriate units.
+
+    Examples
+    --------
+    >>> import astropy.units as u
+    >>> T_e = 5000 * u.K
+    >>> n = 1e19 * u.m ** -3
+    >>> ionization_balance(n, T_e)
+    <Quantity 0.27478417>
+    >>> T_e = 50 * u.eV
+    >>> n = 1e10 * u.m ** -3
+    >>> ionization_balance(n, T_e)
+    <Quantity 12.61576334>
+
+    Returns
+    -------
+    Z: `~astropy.units.Quantity`
+        The average ionization state of the ions in the plasma
+        assuming that the number of ions in each state are equal.
+
+    """
+    E_H = 1 * u.Ry
+
+    A = sqrt(k_B * T_e / E_H)
+    B = log(1 / (4 * n * a0 ** 3) * (k_B * T_e / (pi * E_H)) ** (3 / 2))
+
+    Z = A * sqrt(B) - 1 / 2
+    return Z
+
+
+Z_bal = ionization_balance
+"""alias for :func:`ionization_balance`"""
diff --git a/plasmapy/formulary/parameters.py b/plasmapy/formulary/parameters.py
index 5068e9b9..e1db4115 100644
--- a/plasmapy/formulary/parameters.py
+++ b/plasmapy/formulary/parameters.py
@@ -412,6 +412,16 @@ def ion_sound_speed(
     return V_S
 
 
+# This dictionary defines coefficients for thermal speeds
+# calculated for different methods and values of ndim.
+# Created here to avoid re-instantiating on each call
+_coefficients = {
+    1: {"most_probable": 0, "rms": 1, "mean_magnitude": 2 / np.pi},
+    2: {"most_probable": 1, "rms": 2, "mean_magnitude": np.pi / 2},
+    3: {"most_probable": 2, "rms": 3, "mean_magnitude": 8 / np.pi},
+}
+
+
 @check_relativistic
 @validate_quantities(
     T={"can_be_negative": False, "equivalencies": u.temperature_energy()},
@@ -423,6 +433,7 @@ def thermal_speed(
     particle: particles.Particle = "e-",
     method="most_probable",
     mass: u.kg = np.nan * u.kg,
+    ndim=3,
 ) -> u.m / u.s:
     r"""
     Return the most probable speed for a particle within a Maxwellian
@@ -448,6 +459,10 @@ def thermal_speed(
         anything, but if set, overrides mass acquired from `particle`. Useful
         with relative velocities of particles.
 
+    ndim : int
+        Dimensionality of space in which to calculate thermal velocity. Valid
+        values are 1,2,3.
+
     Returns
     -------
     V : ~astropy.units.Quantity
@@ -479,10 +494,34 @@ def thermal_speed(
     The particle thermal speed is given by:
 
     .. math::
-        V_{th,i} = \sqrt{\frac{2 k_B T_i}{m_i}}
-
-    This function yields the most probable speed within a distribution
-    function.  However, the definition of thermal velocity varies by
+        V_{th,i} = \sqrt{\frac{N k_B T_i}{m_i}}
+
+    where the value of N depends on the dimensionality and the definition of
+    :math:`v_{th}`: most probable, root-mean-square (RMS), or mean magnitude.
+    The value of N in each case is
+
+    .. list-table:: Values of constant N
+       :widths: 50, 25, 25, 25
+       :header-rows: 1
+
+       * - Dim.
+         - Most-Probable
+         - RMS
+         - Mean-Magnitude
+       * - 1D
+         - 0
+         - 1
+         - :math:`2/\pi`
+       * - 2D
+         - 1
+         - 2
+         - :math:`\pi/2`
+       * - 3D
+         - 2
+         - 3
+         - :math:`8/\pi`
+
+    The definition of thermal velocity varies by
     the square root of two depending on whether or not this velocity
     absorbs that factor in the expression for a Maxwellian
     distribution.  In particular, the expression given in the NRL
@@ -509,16 +548,16 @@ def thermal_speed(
     m = mass if np.isfinite(mass) else particles.particle_mass(particle)
 
     # different methods, as per https://en.wikipedia.org/wiki/Thermal_velocity
-    if method == "most_probable":
-        V = np.sqrt(2 * k_B * T / m)
-    elif method == "rms":
-        V = np.sqrt(3 * k_B * T / m)
-    elif method == "mean_magnitude":
-        V = np.sqrt(8 * k_B * T / (m * np.pi))
-    else:
+    try:
+        coef = _coefficients[ndim]
+    except KeyError:
+        raise ValueError("{ndim} is not a supported value for ndim in " "thermal_speed")
+    try:
+        coef = coef[method]
+    except KeyError:
         raise ValueError("Method {method} not supported in thermal_speed")
 
-    return V
+    return np.sqrt(coef * k_B * T / m)
 
 
 @validate_quantities(
diff --git a/requirements/requirements.txt b/requirements/requirements.txt
index 0a968d30..62af843d 100644
--- a/requirements/requirements.txt
+++ b/requirements/requirements.txt
@@ -5,6 +5,6 @@ h5py >= 2.8
 lmfit >= 1.0.1
 matplotlib >= 2.0
 mpmath >= 1.0.0
-numpy >= 1.16.6
+numpy >= 1.18.1
 scipy >= 1.2
 setuptools >= 41.2
diff --git a/setup.cfg b/setup.cfg
index 7f25475c..d8b3f845 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -32,7 +32,7 @@ setup_requires = setuptools_scm
 install_requires =
   astropy >= 4.0
   colorama >= 0.3
-  numpy >= 1.16.6
+  numpy >= 1.18.1
   scipy >= 1.2
   setuptools >= 41.2
 
diff --git a/tox.ini b/tox.ini
index e090249c..c46a57a4 100644
--- a/tox.ini
+++ b/tox.ini
@@ -45,7 +45,7 @@ deps =
   pytest-cov
   pytest-xdist
 conda_deps =
-  numpy>=1.16.6
+  numpy>=1.18.1
   scipy>=1.2
   astropy>=4.0
   colorama>=0.3
@@ -68,7 +68,7 @@ deps =
   pytest-cov
   pytest-xdist
 conda_deps =
-  numpy=1.16.6
+  numpy=1.18.1
   scipy=1.2
   astropy=4.0
   colorama=0.3


Test patch:
diff --git a/plasmapy/diagnostics/tests/test_thomson.py b/plasmapy/diagnostics/tests/test_thomson.py
new file mode 100644
index 00000000..e16b095e
--- /dev/null
+++ b/plasmapy/diagnostics/tests/test_thomson.py
@@ -0,0 +1,251 @@
+"""
+Tests for Thomson scattering analysis functions
+"""
+
+import astropy.units as u
+import numpy as np
+import pytest
+
+from plasmapy.diagnostics import thomson
+from plasmapy.particles import Particle
+
+
+def width_at_value(x, y, val):
+    """
+    Calculates the width of a curve at a given value.
+    """
+    above = np.where(y > val, x, np.nan)
+    return np.abs(np.nanmax(above) - np.nanmin(above))
+
+
+def gen_collective_spectrum():
+    """
+    Generates an example Thomson scattering spectrum in the collective regime
+    """
+    wavelengths = np.arange(520, 545, 0.01) * u.nm
+    probe_wavelength = 532 * u.nm
+    ne = 5e17 * u.cm ** -3
+    probe_vec = np.array([1, 0, 0])
+    scatter_vec = np.array([0, 1, 0])
+    fract = np.array([1.0])
+    Te = 10 * u.eV
+    Ti = np.array([10]) * u.eV
+    ion_species = ["C-12 5+"]
+
+    alpha, Skw = thomson.spectral_density(
+        wavelengths,
+        probe_wavelength,
+        ne,
+        Te,
+        Ti,
+        fract=fract,
+        ion_species=ion_species,
+        probe_vec=probe_vec,
+        scatter_vec=scatter_vec,
+    )
+
+    return alpha, wavelengths, Skw
+
+
+def gen_multiple_ion_species_spectrum():
+    """
+    Generates an example Thomson scattering spectrum for multiple ion species
+    that also have drift velocities. Parameters are set to be in the
+    collective regime where ion species are important.
+    """
+    wavelengths = np.arange(520, 545, 0.01) * u.nm
+    probe_wavelength = 532 * u.nm
+    ne = 5e17 * u.cm ** -3
+    probe_vec = np.array([1, 0, 0])
+    scatter_vec = np.array([0, 1, 0])
+    fract = np.array([0.7, 0.3])
+    Te = 10 * u.eV
+    Ti = np.array([5, 5]) * u.eV
+    fluid_vel = np.array([300, 0, 0]) * u.km / u.s
+    ion_vel = np.array([[-500, 0, 0], [0, 500, 0]]) * u.km / u.s
+
+    # Use this to also test passing in ion species as Particle objects
+    ion_species = [Particle("p+"), Particle("C-12 5+")]
+
+    alpha, Skw = thomson.spectral_density(
+        wavelengths,
+        probe_wavelength,
+        ne,
+        Te,
+        Ti,
+        fract=fract,
+        ion_species=ion_species,
+        probe_vec=probe_vec,
+        scatter_vec=scatter_vec,
+        fluid_vel=fluid_vel,
+        ion_vel=ion_vel,
+    )
+
+    return alpha, wavelengths, Skw
+
+
+def gen_non_collective_spectrum():
+    """
+    Generates an example Thomson scattering spectrum in the non-collective
+    regime
+    """
+    wavelengths = np.arange(500, 570, 0.01) * u.nm
+    probe_wavelength = 532 * u.nm
+    ne = 5e15 * u.cm ** -3
+    probe_vec = np.array([1, 0, 0])
+    scatter_vec = np.array([0, 1, 0])
+    fract = np.array([1.0])
+    Te = 100 * u.eV
+    Ti = np.array([10]) * u.eV
+    ion_species = ["H+"]
+
+    alpha, Skw = thomson.spectral_density(
+        wavelengths,
+        probe_wavelength,
+        ne,
+        Te,
+        Ti,
+        fract=fract,
+        ion_species=ion_species,
+        probe_vec=probe_vec,
+        scatter_vec=scatter_vec,
+    )
+
+    return alpha, wavelengths, Skw
+
+
+def test_different_input_types():
+
+    # Define some constants
+    wavelengths = np.arange(520, 545, 0.01) * u.nm
+    probe_wavelength = 532 * u.nm
+    ne = 5e17 * u.cm ** -3
+    probe_vec = np.array([1, 0, 0])
+    scatter_vec = np.array([0, 1, 0])
+    fract = np.array([1.0])
+    Te = 10 * u.eV
+    Ti = np.array([10]) * u.eV
+    ion_species = "C-12 5+"
+
+    # Raise a ValueError with inconsistent ion array lengths
+    with pytest.raises(ValueError):
+        alpha, Skw = thomson.spectral_density(
+            wavelengths,
+            probe_wavelength,
+            ne,
+            Te,
+            Ti,
+            fract=np.array([0.5, 0.5]),
+            ion_species=ion_species,
+            probe_vec=probe_vec,
+            scatter_vec=scatter_vec,
+        )
+
+    # Raise a ValueError with inconsistent ion temperature array
+    with pytest.raises(ValueError):
+        alpha, Skw = thomson.spectral_density(
+            wavelengths,
+            probe_wavelength,
+            ne,
+            Te,
+            np.array([5, 5]) * u.eV,
+            fract=fract,
+            ion_species=ion_species,
+            probe_vec=probe_vec,
+            scatter_vec=scatter_vec,
+        )
+
+    # Raise a ValueError with empty ion_species
+    with pytest.raises(ValueError):
+        alpha, Skw = thomson.spectral_density(
+            wavelengths,
+            probe_wavelength,
+            ne,
+            Te,
+            Ti,
+            fract=fract,
+            ion_species=[],
+            probe_vec=probe_vec,
+            scatter_vec=scatter_vec,
+        )
+
+
+def test_collective_spectrum():
+    """
+    Compares the generated spectrum to previously determined values
+    """
+    alpha, wavelength, Skw = gen_collective_spectrum()
+
+    # Check that alpha is correct
+    assert np.isclose(alpha.value, 1.801, atol=0.01), (
+        "Collective case alpha " f"returns {alpha} instead of " "expected 1.801"
+    )
+
+    i_width = width_at_value(wavelength.value, Skw.value, 2e-13)
+    e_width = width_at_value(wavelength.value, Skw.value, 0.2e-13)
+
+    # Check that the widths of the ion and electron features match expectations
+    assert np.isclose(i_width, 0.1599, 1e-3), (
+        "Collective case ion feature "
+        f"width is {i_width}"
+        "instead of expected 0.1599"
+    )
+
+    assert np.isclose(e_width, 17.7899, 1e-3), (
+        "Collective case electron "
+        f"feature width is {e_width} "
+        "instead of expected 17.7899"
+    )
+
+
+def test_non_collective_spectrum():
+    """
+    Compares the generated spectrum to previously determined values
+    """
+    alpha, wavelength, Skw = gen_non_collective_spectrum()
+
+    # Check that alpha is correct
+    assert np.isclose(alpha.value, 0.05707, atol=0.01), (
+        "Non-collective case alpha " f"returns {alpha} instead of " "expected 0.05707"
+    )
+
+    e_width = width_at_value(wavelength.value, Skw.value, 0.2e-13)
+
+    # Check that the widts of the electron feature matchs expectations
+    assert np.isclose(e_width, 22.6699, 1e-3), (
+        "Non-collective case electron "
+        f"feature width is {e_width} "
+        "instead of expected 22.6699"
+    )
+
+
+def test_multiple_ion_species_spectrum():
+    """
+    Compares the generated spectrum to previously determined values
+    """
+
+    alpha, wavelength, Skw = gen_multiple_ion_species_spectrum()
+
+    # Compute the width and max of the spectrum, and the wavelength
+    # of the max (sensitive to ion vel)
+    width = width_at_value(wavelength.value, Skw.value, 0.2e-11)
+    max_skw = np.max(Skw.value)
+    max_wavelength = wavelength.value[np.argmax(Skw.value)]
+
+    # Check width
+    assert np.isclose(width, 0.1599, 1e-3), (
+        f"Multiple ion species case spectrum width is {width} instead of "
+        "expected 0.1599"
+    )
+
+    # Check max value
+    assert np.isclose(max_skw, 2.4e-11, 1e-11), (
+        f"Multiple ion species case spectrum max is {max_skw} instead of "
+        "expected 2.4e-11"
+    )
+
+    # Check max peak location
+    assert np.isclose(max_wavelength, 531.549, 1e-3), (
+        "Multiple ion species case spectrum peak wavelength is "
+        f"{max_wavelength} instead of expected 531.549"
+    )
diff --git a/plasmapy/formulary/tests/test_ionization.py b/plasmapy/formulary/tests/test_ionization.py
new file mode 100644
index 00000000..3df0aa45
--- /dev/null
+++ b/plasmapy/formulary/tests/test_ionization.py
@@ -0,0 +1,19 @@
+import astropy.units as u
+import numpy as np
+import pytest
+from astropy.constants import c
+
+from plasmapy.formulary.ionization import Z_bal
+
+
+def test_Z_bal():
+    n = 1e19 * u.m ** -3
+    T_e = 5000 * u.K
+
+    assert (Z_bal(n, T_e) * u.dimensionless_unscaled).unit == u.dimensionless_unscaled
+
+    with pytest.warns(u.UnitsWarning):
+        Z_bal(1e19, T_e)
+
+    with pytest.raises(u.UnitTypeError):
+        Z_bal(2e19 * u.kg, T_e)
diff --git a/plasmapy/formulary/tests/test_parameters.py b/plasmapy/formulary/tests/test_parameters.py
index 7fc8486b..bc5435f8 100644
--- a/plasmapy/formulary/tests/test_parameters.py
+++ b/plasmapy/formulary/tests/test_parameters.py
@@ -365,6 +365,45 @@ def test_thermal_speed():
         thermal_speed(1 * u.MK, particle="p").si.value, 128486.56960876315
     )
 
+    # Explicitly check all three modes and dimensionalities
+    # ndim = 1
+    assert np.isclose(thermal_speed(T_e, method="most_probable", ndim=1).si.value, 0.0)
+
+    # Regression tests start here!
+    assert np.isclose(
+        thermal_speed(T_e, method="rms", ndim=1).si.value, 3893114.2008620175
+    )
+
+    assert np.isclose(
+        thermal_speed(T_e, method="mean_magnitude", ndim=1).si.value, 3106255.714310189
+    )
+
+    # ndim = 2
+    assert np.isclose(
+        thermal_speed(T_e, method="most_probable", ndim=2).si.value, 3893114.2008620175
+    )
+
+    assert np.isclose(
+        thermal_speed(T_e, method="rms", ndim=2).si.value, 5505694.902726359
+    )
+
+    assert np.isclose(
+        thermal_speed(T_e, method="mean_magnitude", ndim=2).si.value, 4879295.066124102
+    )
+
+    # ndim = 3
+    assert np.isclose(
+        thermal_speed(T_e, method="most_probable", ndim=3).si.value, 5505694.902726359
+    )
+
+    assert np.isclose(
+        thermal_speed(T_e, method="rms", ndim=3).si.value, 6743071.595560921
+    )
+
+    assert np.isclose(
+        thermal_speed(T_e, method="mean_magnitude", ndim=3).si.value, 6212511.428620378
+    )
+
     # Case when Z=1 is assumed
     assert thermal_speed(T_i, particle="p") == thermal_speed(T_i, particle="H-1+")
 
@@ -398,9 +437,14 @@ def test_thermal_speed():
         thermal_speed(1e6 * u.K, method="rms").si.value, 6743070.475775486
     )
 
+    # Test invalid method
     with pytest.raises(ValueError):
         thermal_speed(T_i, method="sadks")
 
+    # Test invalid ndim
+    with pytest.raises(ValueError):
+        thermal_speed(T_i, ndim=4)
+
     assert_can_handle_nparray(thermal_speed)