Repository: PlasmaPy/PlasmaPy
Pull request #893 (2020-08-26 17:57:14)

Thomson Scattering Multiple Electron Populations (and a bug fix)
This PR adds support to the diagnostics.Thomson.spectral_density() function for multiple thermal electron populations (still under the requirement that the resulting plasma is quasi-neutral). The electron keywords are adjusted to use the same syntax as the ion keywords (the function already supported multiple ion species).

This PR also fixes a bug for multiple ion species by adding a factor of ifract to the calculation of the ion contribution.  Physically, this corresponds to modifying the magnitude of the Maxwellian distribution function in S(k,w). Without this, running the code for two identical ion populations with ifract=[0.5,0.5] gives twice the result as running it for ifract=[1.0]. The same correction factor has been added for efract and both now pass this test (which has been added explicitly as a unit test).

review-notebook-app[bot]: Check out this pull request on&nbsp; <a href="https://app.reviewnb.com/PlasmaPy/PlasmaPy/pull/893"><img align="absmiddle"  alt="ReviewNB" height="28" class="BotMessageButtonImage" src="https://raw.githubusercontent.com/ReviewNB/support/master/images/button_reviewnb.png"/></a> 

 Review Jupyter notebook visual diffs & provide feedback on notebooks. 

---

 <i>Powered by <a href='https://www.reviewnb.com/?utm_source=gh'>ReviewNB</a></i>
codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/893?src=pr&el=h1) Report
> Merging [#893](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/893?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/37120c7cf0c0703b3363e17d678746c6d9ed6f76?el=desc) will **decrease** coverage by `0.06%`.
> The diff coverage is `93.33%`.

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

```diff
@@            Coverage Diff             @@
##           master     #893      +/-   ##
==========================================
- Coverage   96.18%   96.11%   -0.07%     
==========================================
  Files          59       59              
  Lines        5367     5333      -34     
==========================================
- Hits         5162     5126      -36     
- Misses        205      207       +2     
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/893?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/diagnostics/thomson.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/893/diff?src=pr&el=tree#diff-cGxhc21hcHkvZGlhZ25vc3RpY3MvdGhvbXNvbi5weQ==) | `97.43% <93.33%> (-2.57%)` | :arrow_down: |
| [plasmapy/plasma/plasma\_base.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/893/diff?src=pr&el=tree#diff-cGxhc21hcHkvcGxhc21hL3BsYXNtYV9iYXNlLnB5) | `72.22% <0.00%> (-3.39%)` | :arrow_down: |
| [plasmapy/simulation/abstractions.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/893/diff?src=pr&el=tree#diff-cGxhc21hcHkvc2ltdWxhdGlvbi9hYnN0cmFjdGlvbnMucHk=) | `75.00% <0.00%> (-2.78%)` | :arrow_down: |
| [plasmapy/plasma/sources/openpmd\_hdf5.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/893/diff?src=pr&el=tree#diff-cGxhc21hcHkvcGxhc21hL3NvdXJjZXMvb3BlbnBtZF9oZGY1LnB5) | `78.88% <0.00%> (-1.12%)` | :arrow_down: |
| [plasmapy/simulation/particletracker.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/893/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/893/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/893/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/893/diff?src=pr&el=tree#diff-cGxhc21hcHkvcGFydGljbGVzL3BhcnRpY2xlX2NsYXNzLnB5) | `96.57% <0.00%> (-0.01%)` | :arrow_down: |
| [plasmapy/formulary/braginskii.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/893/diff?src=pr&el=tree#diff-cGxhc21hcHkvZm9ybXVsYXJ5L2JyYWdpbnNraWkucHk=) | `99.72% <0.00%> (-0.01%)` | :arrow_down: |
| [plasmapy/particles/decorators.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/893/diff?src=pr&el=tree#diff-cGxhc21hcHkvcGFydGljbGVzL2RlY29yYXRvcnMucHk=) | `100.00% <0.00%> (ø)` | |
| ... and [6 more](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/893/diff?src=pr&el=tree-more) | |

------

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

pheuer: @StanczakDominik Thank you for the edits! I've made a slight change to the notebook example to include multiple electron populations in the last set of plots (along with multiple ion populations).

Patch:
diff --git a/changelog/892.feature.rst b/changelog/892.feature.rst
new file mode 100644
index 00000000..22cea5dc
--- /dev/null
+++ b/changelog/892.feature.rst
@@ -0,0 +1,1 @@
+Added a function to calculate the power spectrum of thermal bremsstrahlung emitted by a Maxwellian plasma.
diff --git a/changelog/893.feature.rst b/changelog/893.feature.rst
new file mode 100644
index 00000000..b7d4406c
--- /dev/null
+++ b/changelog/893.feature.rst
@@ -0,0 +1,1 @@
+Added support for multiple electron components to diagnostics.thomson.spectral_density. Also fixed a bug for multiple ion populations.
diff --git a/docs/formulary/index.rst b/docs/formulary/index.rst
index dbf3cd97..90f3f9c2 100644
--- a/docs/formulary/index.rst
+++ b/docs/formulary/index.rst
@@ -49,6 +49,9 @@ physical quantities helpful for plasma physics.
    | .. toctree:: Quantum Relations <quantum>               | `plasmapy.formulary.quantum`            |
    |    :maxdepth: 1                                        |                                         |
    +--------------------------------------------------------+-----------------------------------------+
+   | .. toctree:: Radiation <radiation>                     | `plasmapy.formulary.radiation`          |
+   |    :maxdepth: 1                                        |                                         |
+   +--------------------------------------------------------+-----------------------------------------+
    | .. toctree:: Relativistic Relations <relativity>       | `plasmapy.formulary.relativity`         |
    |    :maxdepth: 1                                        |                                         |
    +--------------------------------------------------------+-----------------------------------------+
diff --git a/docs/formulary/radiation.rst b/docs/formulary/radiation.rst
new file mode 100644
index 00000000..6272ce2a
--- /dev/null
+++ b/docs/formulary/radiation.rst
@@ -0,0 +1,16 @@
+.. _radiation:
+
+************************************************************************
+Electromagnetic Radiation Functions (`plasmapy.formulary.radiation`)
+************************************************************************
+
+.. currentmodule:: plasmapy.formulary.radiation
+
+.. automodapi:: plasmapy.formulary.radiation
+   :no-heading:
+
+
+.. nbgallery::
+    :caption: Examples
+
+    /notebooks/thermal_bremsstrahlung
diff --git a/docs/notebooks/thermal_bremsstrahlung.ipynb b/docs/notebooks/thermal_bremsstrahlung.ipynb
new file mode 100644
index 00000000..08320c7a
--- /dev/null
+++ b/docs/notebooks/thermal_bremsstrahlung.ipynb
@@ -0,0 +1,170 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Emission of Thermal Bremsstrahlung by a Maxwellian Plasma"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The radiation.thermal_bremsstrahlung function calculates the bremsstrahlung spectrum emitted by the collision of electrons and ions in a thermal (Maxwellian) plasma. This function calculates this quantity in the Rayleigh-Jeans limit where $\\hbar\\omega \\ll k_B T_e$. In this regime, the power spectrum of the emitted radiation is\n",
+    "\n",
+    "\\begin{equation}\n",
+    "\\frac{dP}{d\\omega} = \\frac{8 \\sqrt{2}}{3\\sqrt{\\pi}} \\bigg ( \\frac{e^2}{4 \\pi \\epsilon_0} \\bigg )^3 \\bigg ( m_e c^2 \\bigg )^{-\\frac{3}{2}} \\bigg ( 1 - \\frac{\\omega_{pe}^2}{\\omega^2} \\bigg )^\\frac{1}{2} \\frac{Z_i^2 n_i n_e}{\\sqrt{k_B T_e}} E_1(y)\n",
+    "\\end{equation}\n",
+    "\n",
+    "where $w_{pe}$ is the electron plasma frequency and $E_1$ is the exponential integral\n",
+    "\n",
+    "\\begin{equation}\n",
+    "E_1 (y) = - \\int_{-y}^\\infty \\frac{e^{-t}}{t}dt\n",
+    "\\end{equation}\n",
+    "\n",
+    "and y is the dimensionless argument\n",
+    "\n",
+    "\\begin{equation}\n",
+    "y = \\frac{1}{2} \\frac{\\omega^2 m_e}{k_{max}^2 k_B T_e}\n",
+    "\\end{equation}\n",
+    "\n",
+    "where $k_{max}$ is a maximum wavenumber arising from binary collisions approximated here as\n",
+    "\n",
+    "\\begin{equation}\n",
+    "k_{max} = \\frac{1}{\\lambda_B} = \\frac{\\sqrt{m_e k_B T_e}}{\\hbar}\n",
+    "\\end{equation}\n",
+    "\n",
+    "where $\\lambda_B$ is the electron de Broglie wavelength. In some regimes other values for $k_{max}$ may be appropriate, so its value may be set using a keyword. Bremsstrahlung emission is greatly reduced below the electron plasma frequency (where the plasma is opaque to EM radiation), so these expressions are only valid in the regime $w < w_{pe}$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "\n",
+    "import astropy.constants as const\n",
+    "import astropy.units as u\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "from plasmapy.formulary.radiation import thermal_bremsstrahlung"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Create an array of frequencies over which to calculate the bremsstrahlung spectrum and convert these frequencies to photon energies for the purpose of plotting the results. Set the plasma density, temperature, and ion species."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "frequencies = np.arange(15, 16, .01)\n",
+    "frequencies = (10**frequencies)/u.s\n",
+    "\n",
+    "energies = (frequencies*const.h.si).to(u.eV)\n",
+    "\n",
+    "ne = 1e22*u.cm**-3\n",
+    "Te = 1e2*u.eV\n",
+    "ion_species = 'C-12 4+'"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Calculate the spectrum, then plot it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "nbsphinx-thumbnail": {
+     "tooltip": "Thermal Bremsstrahlung"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "spectrum = thermal_bremsstrahlung(frequencies, ne, Te, ion_species=ion_species)\n",
+    "\n",
+    "print(spectrum.unit)\n",
+    "\n",
+    "lbl = \"$T_e$ = {:.1e} eV,\\n\".format(Te.value) + \"$n_e$ = {:.1e} 1/cm^3\".format(ne.value)\n",
+    "plt.plot(energies, spectrum, label=lbl)\n",
+    "plt.title(f\"Thermal Bremsstrahlung Spectrum\")\n",
+    "plt.xlabel('Energy (eV)')\n",
+    "plt.ylabel('Power Spectral Density (W s/m^3)')\n",
+    "plt.legend()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The power spectrum is the power per angular frequency per volume integrated over $4\\pi$ sr of solid angle, and therefore has units of Watts / (rad/s) / m$^3$ * $4\\pi$ rad = W s/m$^3$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "spectrum = spectrum.to(u.W*u.s/u.m**3)\n",
+    "spectrum.unit"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This means that, for a given volume and time period, the total energy emitted can be determined by integrating the power spectrum"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "t = 5 * u.ns\n",
+    "vol = 0.5*u.cm**3\n",
+    "dw = 2*np.pi*np.gradient(frequencies) # Frequency step size\n",
+    "total_energy =  (np.sum(spectrum*dw)*t*vol).to(u.J)\n",
+    "print(\"Total Energy: {:.2e} J\".format(total_energy.value))"
+   ]
+  }
+ ],
+ "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
index 70844290..624430d5 100644
--- a/docs/notebooks/thomson.ipynb
+++ b/docs/notebooks/thomson.ipynb
@@ -15,13 +15,13 @@
     "[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",
+    "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 electron populations under the assumption that all of the ion species and the electron fluid have Maxwellian velocity distribution functions and that the combined plasma is quasi-neutral. 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",
+    "S(k,\\omega) = \\sum_e \\frac{2\\pi}{k} \\bigg |1 - \\frac{\\chi_e}{\\epsilon} \\bigg |^2 f_{e0,e}\\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",
+    "where $\\chi_e$ is the electron component susceptibility of the plasma and $\\epsilon = 1 + \\sum_e \\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,e}$ 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",
@@ -135,12 +135,11 @@
    "source": [
     "non_collective = {\n",
     "    'name': 'Non-Collective Regime', \n",
-    "    'ne': 5e15*u.cm**-3, \n",
-    "    'Te': 40*u.eV, \n",
+    "    'n': 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",
+    "    'electron_vel': np.array([[0, 0, 0]])*u.km/u.s, \n",
     "    'ion_vel': np.array([[0, 0, 0]])*u.km/u.s,\n",
     "}"
    ]
@@ -160,12 +159,11 @@
    "source": [
     "weakly_collective = {\n",
     "    'name': 'Weakly Collective Regime',\n",
-    "    'ne': 2e17*u.cm**-3,\n",
+    "    'n': 2e17*u.cm**-3,\n",
     "    'Te': 20*u.eV,\n",
-    "    'Ti': np.array([10])*u.eV,\n",
-    "    'fract': np.array([1]),\n",
+    "    'Ti': 10*u.eV,\n",
     "    'ion_species': ['H+'],\n",
-    "    'fluid_vel': np.array([0, 0, 0])*u.km/u.s,\n",
+    "    'electron_vel': np.array([[0, 0, 0]])*u.km/u.s,\n",
     "    'ion_vel': np.array([[0, 0, 0]])*u.km/u.s,\n",
     "}"
    ]
@@ -185,12 +183,11 @@
    "source": [
     "collective = {\n",
     "    'name': 'Collective Regime',\n",
-    "    'ne': 5e17*u.cm**-3,\n",
+    "    'n': 5e17*u.cm**-3,\n",
     "    'Te': 10*u.eV,\n",
-    "    'Ti': np.array([4])*u.eV,\n",
-    "    'fract': np.array([1]),\n",
+    "    'Ti': 4*u.eV,\n",
     "    'ion_species': ['H+'],\n",
-    "    'fluid_vel': np.array([0, 0, 0])*u.km/u.s,\n",
+    "    'electron_vel': np.array([[0, 0, 0]])*u.km/u.s,\n",
     "    'ion_vel': np.array([[0, 0, 0]])*u.km/u.s,\n",
     "}"
    ]
@@ -210,12 +207,11 @@
    "source": [
     "drifts = {\n",
     "    'name': 'Drift Velocities',\n",
-    "    'ne': 5e17*u.cm**-3,\n",
+    "    'n': 5e17*u.cm**-3,\n",
     "    'Te': 10*u.eV,\n",
-    "    'Ti': np.array([10])*u.eV,\n",
-    "    'fract': np.array([1]),\n",
+    "    'Ti': 10*u.eV,\n",
     "    'ion_species': ['H+'],\n",
-    "    'fluid_vel': np.array([700, 0, 0])*u.km/u.s,\n",
+    "    'electron_vel': np.array([[700, 0, 0]])*u.km/u.s,\n",
     "    'ion_vel': np.array([[-600, -100, 0]])*u.km/u.s,\n",
     "}"
    ]
@@ -224,7 +220,7 @@
    "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"
+    "A case identical to the collective example, except that now the plasma consists 25% He+1 and 75% C+5, and two electron populations exist with different temperatures."
    ]
   },
   {
@@ -234,13 +230,14 @@
    "outputs": [],
    "source": [
     "two_species = {\n",
-    "    'name': 'Two Ion Species',\n",
-    "    'ne': 5e17*u.cm**-3,\n",
-    "    'Te': 10*u.eV,\n",
+    "    'name': 'Two Ion and Electron Components',\n",
+    "    'n': 5e17*u.cm**-3,\n",
+    "    'Te': np.array([50, 10])*u.eV,\n",
     "    'Ti': np.array([10, 50])*u.eV,\n",
-    "    'fract': np.array([.25, .75]),\n",
+    "    'efract': np.array([0.5, 0.5]), \n",
+    "    'ifract': 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",
+    "    'electron_vel': np.array([[0, 0, 0],[0, 0, 0]])*u.km/u.s,\n",
     "    'ion_vel': np.array([[0, 0, 0],[0, 0, 0]])*u.km/u.s,\n",
     "}"
    ]
@@ -275,8 +272,8 @@
     "for i, x in enumerate(examples):\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",
+    "                     x['n'], x['Te'], x['Ti'], ifract=x.get('ifract'), efract=x.get('efract'),\n",
+    "                     ion_species=x['ion_species'], electron_vel=x['electron_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",
@@ -310,8 +307,15 @@
     "\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"
+    "e. Including multiple ion and electron populations modifies the ion and electron features respectively.\n"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
diff --git a/plasmapy/diagnostics/thomson.py b/plasmapy/diagnostics/thomson.py
index 3a4d95c5..aa6c82a4 100644
--- a/plasmapy/diagnostics/thomson.py
+++ b/plasmapy/diagnostics/thomson.py
@@ -27,19 +27,20 @@
 @validate_quantities(
     wavelengths={"can_be_negative": False},
     probe_wavelength={"can_be_negative": False},
-    ne={"can_be_negative": False},
+    n={"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,
+    n: u.m ** -3,
     Te: u.K,
     Ti: u.K,
-    fract: np.ndarray = np.ones(1),
+    efract: np.ndarray = None,
+    ifract: np.ndarray = None,
     ion_species: Union[str, List[str], Particle, List[Particle]] = "H+",
-    fluid_vel: u.m / u.s = np.zeros(3) * u.m / u.s,
+    electron_vel: u.m / u.s = None,
     ion_vel: u.m / u.s = None,
     probe_vec=np.array([1, 0, 0]),
     scatter_vec=np.array([0, 1, 0]),
@@ -50,22 +51,23 @@ def spectral_density(
 
     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:
+    species and a one or more thermal electron populations (the entire plasma
+    is assumed to be quasi-neutral)
 
     .. math::
-        S(k,\omega) = 1 - \frac{2\pi}{k}
+        S(k,\omega) = \sum_e \frac{2\pi}{k}
         \bigg |1 - \frac{\chi_e}{\epsilon} \bigg |^2
-        f_{e0} \bigg (\frac{\omega}{k} \bigg ) +
+        f_{e0,e} \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 and :math:`\epsilon = 1 + \sum_e \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
+    and :math:`f_{e0,e}` 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`_.
@@ -80,34 +82,40 @@ def spectral_density(
     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.
+    n : `~astropy.units.Quantity`
+        Mean (0th order) 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)
+    Te : `~astropy.units.Quantity`, shape (Ne, )
+        Temperature of each electron component. Shape (Ne, ) must be equal to the
+        number of electron components Ne. (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)
+    Ti : `~astropy.units.Quantity`, shape (Ni, )
+        Temperature of each ion component. Shape (Ni, ) must be equal to the
+        number of ion components Ni. (in K or convertible to eV)
 
-    fract : array_like, shape (N, ), optional
+    efract : array_like, shape (Ne, ), optional
+        An array-like object where each element represents the fraction (or ratio)
+        of the electron component number density to the total electron number density.
+        Must sum to 1.0. Default is a single electron component.
+
+    ifract : array_like, shape (Ni, ), 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
+    ion_species : str or `~plasmapy.particles.Particle`, shape (Ni, ), 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)
+    electron_vel : `~astropy.units.Quantity`, shape (Ne, 3), optional
+        Velocity of each electron component 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
+    ion_vel : `~astropy.units.Quantity`, shape (Ni, 3), optional
+        Velocity vectors for each electron population in the rest frame
+        (convertible to m/s) Defaults zero drift
         for all specified ion species.
 
     probe_vec : float `~numpy.ndarray`, shape (3, )
@@ -123,7 +131,8 @@ def spectral_density(
     -------
     alpha : float
         Mean scattering parameter, where `alpha` > 1 corresponds to collective
-        scattering and `alpha` < 1 indicates non-collective scattering.
+        scattering and `alpha` < 1 indicates non-collective scattering. The
+        scattering parameter is calculated based on the total plasma density n.
 
     Skw : `~astropy.units.Quantity`
         Computed spectral density function over the input `wavelengths` array
@@ -143,12 +152,25 @@ def spectral_density(
     .. _`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 efract is None:
+        efract = np.ones(1)
+    else:
+        efract = np.asarray(efract, dtype=np.float64)
+
+    if ifract is None:
+        ifract = np.ones(1)
+    else:
+        ifract = np.asarray(ifract, dtype=np.float64)
+
+    # If electon velocity is not specified, create an array corresponding
+    # to zero drift
+    if electron_vel is None:
+        electron_vel = np.zeros([efract.size, 3]) * u.m / u.s
 
     # 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
+        ion_vel = np.zeros([ifract.size, 3]) * u.m / u.s
 
     # Condition ion_species
     if isinstance(ion_species, (str, Particle)):
@@ -160,6 +182,16 @@ def spectral_density(
             continue
         ion_species[ii] = Particle(ion)
 
+    # Condition Te
+    if Te.size == 1:
+        # If a single quantity is given, put it in an array so it's iterable
+        # If Te.size != len(efract), assume same temp. for all species
+        Te = np.repeat(Te, len(efract))
+    elif Te.size != len(efract):
+        raise ValueError(
+            f"Got {Te.size} electron temperatures and expected " f"{len(efract)}."
+        )
+
     # Condition Ti
     if Ti.size == 1:
         # If a single quantity is given, put it in an array so it's iterable
@@ -170,14 +202,23 @@ def spectral_density(
             f"Got {Ti.size} ion temperatures and expected " f"{len(ion_species)}."
         )
 
+    # Make sure the sizes of ion_species, ifract, ion_vel, and Ti all match
     if (
-        (len(ion_species) != fract.size)
-        or (ion_vel.shape[0] != fract.size)
-        or (Ti.size != fract.size)
+        (len(ion_species) != ifract.size)
+        or (ion_vel.shape[0] != ifract.size)
+        or (Ti.size != ifract.size)
     ):
         raise ValueError(
-            "Inconsistent number of species in fract, "
-            "ion_species, Ti, and/or ion_vel."
+            f"Inconsistent number of species in ifract ({ifract}), "
+            f"ion_species ({len(ion_species)}), Ti ({Ti.size}), "
+            f"and/or ion_vel ({ion_vel.shape[0]})."
+        )
+
+    # Make sure the sizes of efract, electron_vel, and Te all match
+    if (electron_vel.shape[0] != efract.size) or (Te.size != efract.size):
+        raise ValueError(
+            f"Inconsistent number of electron populations in efract ({efract.size}), "
+            f"Te ({Te.size}), or electron velocity ({electron_vel.shape[0]})."
         )
 
     # Ensure unit vectors are normalized
@@ -194,12 +235,11 @@ def spectral_density(
         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
+    zbar = np.sum(ifract * ion_z)
+    ne = efract * n
+    ni = ifract * n / zbar  # ne/zbar = sum(ni)
+    # wpe is calculated for the entire plasma (all electron populations combined)
+    wpe = plasma_frequency(n=n, particle="e-")
 
     # Convert wavelengths to angular frequencies (electromagnetic waves, so
     # phase speed is c)
@@ -224,54 +264,58 @@ def spectral_density(
     # 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_e = w - np.matmul(electron_vel, np.outer(k, k_vec).T)
     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)
+    alpha = np.sqrt(2) * wpe / np.outer(k, vTe)
 
     # Calculate the normalized phase velocities (Sec. 3.4.2 in Sheffield)
-    xe = (w_e / (k * vTe)).to(u.dimensionless_unscaled)
+    xe = (np.outer(1 / vTe, 1 / k) * w_e).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-")
+
+    chiE = np.zeros([efract.size, w.size], dtype=np.complex128)
+    for i, fract in enumerate(efract):
+        chiE[i, :] = permittivity_1D_Maxwellian(w_e[i, :], k, Te[i], ne[i], "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)
+    chiI = np.zeros([ifract.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, :] = (
+    epsilon = 1 + np.sum(chiE, axis=0) + np.sum(chiI, axis=0)
+
+    econtr = np.zeros([efract.size, w.size], dtype=np.complex128) * u.s / u.rad
+    for m in range(efract.size):
+        econtr[m, :] = efract[m] * (
+            2
+            * np.sqrt(np.pi)
+            / k
+            / vTe[m]
+            * np.power(np.abs(1 - np.sum(chiE, axis=0) / epsilon), 2)
+            * np.exp(-xe[m, :] ** 2)
+        )
+
+    icontr = np.zeros([ifract.size, w.size], dtype=np.complex128) * u.s / u.rad
+    for m in range(ifract.size):
+        icontr[m, :] = ifract[m] * (
             2
             * np.sqrt(np.pi)
             * ion_z[m]
             / k
             / vTi[m]
-            * np.power(np.abs(chiE / epsilon), 2)
+            * np.power(np.abs(np.sum(chiE, axis=0) / epsilon), 2)
             * np.exp(-xi[m, :] ** 2)
         )
 
     # Recast as real: imaginary part is already zero
-    Skw = np.real(econtr + np.sum(icontr, axis=0))
+    Skw = np.real(np.sum(econtr, axis=0) + np.sum(icontr, axis=0))
 
     return np.mean(alpha), Skw
diff --git a/plasmapy/formulary/radiation.py b/plasmapy/formulary/radiation.py
new file mode 100644
index 00000000..62caa436
--- /dev/null
+++ b/plasmapy/formulary/radiation.py
@@ -0,0 +1,160 @@
+"""
+Defin_es the Radiation module as part of the Plasma Formulary
+"""
+
+__all__ = [
+    "thermal_bremsstrahlung",
+]
+
+import astropy.constants as const
+import astropy.units as u
+import numpy as np
+
+from scipy.special import exp1
+from typing import List, Tuple, Union
+
+from plasmapy.formulary.parameters import plasma_frequency
+from plasmapy.particles import Particle, particle_input
+from plasmapy.utils.decorators import validate_quantities
+from plasmapy.utils.exceptions import PhysicsError
+
+
+@validate_quantities(
+    frequencies={"can_be_negative": False},
+    n_e={"can_be_negative": False},
+    n_i={"can_be_negative": False},
+    T_e={"can_be_negative": False, "equivalencies": u.temperature_energy()},
+)
+@particle_input
+def thermal_bremsstrahlung(
+    frequencies: u.Hz,
+    n_e: u.m ** -3,
+    T_e: u.K,
+    n_i: u.m ** -3 = None,
+    ion_species: Particle = "H+",
+    kmax: u.m = None,
+) -> np.ndarray:
+    r"""
+    Calculate the Bremsstrahlung emission spectrum for a Maxwellian plasma
+    in the Rayleigh-Jeans limit :math:`\hbar\omega \ll k_B*T_e`
+
+    .. math::
+       \frac{dP}{d\omega} = \frac{8 \sqrt{2}}{3\sqrt{\pi}}
+       \bigg ( \frac{e^2}{4 \pi \epsilon_0} \bigg )^3
+       \bigg ( m_e c^2 \bigg )^{-\frac{3}{2}}
+       \bigg ( 1 - \frac{\omega_{pe}^2}{\omega^2} \bigg )^\frac{1}{2}
+       \frac{Z_i^2 n_i n_e}{\sqrt(k_B T_e)}
+       E_1(y)
+
+    where E_1 is the exponential integral
+
+    .. math::
+        E_1 (y) = - \int_{-y}^\infty \frac{e^{-t}}{t}dt
+
+    and y is the dimensionless argument
+
+    .. math::
+        y = \frac{1}{2} \frac{\omega^2 m_e}{k_{max}^2 k_B T_e}
+
+    where   :math:`k_{max}` is a maximum wavenumber approximated here as
+    :math:`k_{max} = 1/\lambda_B` where  :math:`\lambda_B` is the electron
+    de Broglie wavelength.
+
+    Parameters
+    ----------
+
+    frequencies : `~astropy.units.Quantity`
+        Array of frequencies over which the bremsstrahlung spectrum will be
+        calculated (convertable to Hz).
+
+    n_e : `~astropy.units.Quantity`
+        Electron number density in the plasma (convertable to m^-3)
+
+    T_e : `~astropy.units.Quantity`
+        Temperature of the electrons (in K or convertible to eV)
+
+    n_i : `~astropy.units.Quantity` (optional)
+        Ion number density in the plasma (convertable to m^-3). Defaults
+        to the quasi-neutral conditon n_i=n_e/Z.
+
+    ion : str or `~plasmapy.particles.Particle`, optional
+        An instance of `~plasmapy.particles.Particle`, or a string
+        convertible to `~plasmapy.particles.Particle`.
+
+    kmax :  `~astropy.units.Quantity`
+        Cutoff wavenumber (convertable to u.rad/u.m). Defaults to the inverse
+        of the electron de Broglie wavelength.
+
+    Returns
+    -------
+
+    spectrum : `~astropy.units.Quantity`
+        Computed bremsstrahlung spectrum over the frequencies provided.
+
+    Notes
+    -----
+
+    For details, see "Radiation Processes in Plasmas" by
+    Bekefi. `ISBN 978\\-0471063506`_.
+
+    .. _`ISBN 978\\-0471063506`: https://ui.adsabs.harvard.edu/abs/1966rpp..book.....B/abstract
+    """
+
+    # Default n_i is n_e/Z:
+    if n_i is None:
+        n_i = n_e / ion_species.integer_charge
+
+    # Default value of kmax is the electrom thermal de Broglie wavelength
+    if kmax is None:
+        kmax = (np.sqrt(const.m_e.si * const.k_B.si * T_e) / const.hbar.si).to(1 / u.m)
+
+    # Convert frequencies to angular frequencies
+    ω = (frequencies * 2 * np.pi * u.rad).to(u.rad / u.s)
+
+    # Calculate the electron plasma frequency
+    ω_pe = plasma_frequency(n=n_e, particle="e-")
+
+    # Check that all ω < wpe (this formula is only valid in this limit)
+    if np.min(ω) < ω_pe:
+        raise PhysicsError(
+            "Lowest frequency must be larger than the electron "
+            f"plasma frequency {ω_pe:.1e}, but min(ω) = {np.min(ω):.1e}"
+        )
+
+    # Check that the parameters given fall within the Rayleigh-Jeans limit
+    # hω << kT_e
+    rj_const = (
+        np.max(ω) * const.hbar.si / (2 * np.pi * u.rad * const.k_B.si * T_e)
+    ).to(u.dimensionless_unscaled)
+    if rj_const.value > 0.1:
+
+        raise PhysicsError(
+            "Rayleigh-Jeans limit not satisfied: "
+            "hbar*ω/kT_e = {rj_const.value:.2e} > 0.1. "
+            "Try lower ω or higher T_e."
+        )
+
+    # Calculate the bremsstrahlung power spectral density in several steps
+    c1 = (
+        (8 / 3)
+        * np.sqrt(2 / np.pi)
+        * (const.e.si ** 2 / (4 * np.pi * const.eps0.si)) ** 3
+        * 1
+        / (const.m_e.si * const.c.si ** 2) ** 1.5
+    )
+
+    Zi = ion_species.integer_charge
+    c2 = (
+        np.sqrt(1 - ω_pe ** 2 / ω ** 2)
+        * Zi ** 2
+        * n_i
+        * n_e
+        / np.sqrt(const.k_B.si * T_e)
+    )
+
+    # Dimensionless argument for exponential integral
+    arg = 0.5 * ω ** 2 * const.m_e.si / (kmax ** 2 * const.k_B.si * T_e) / u.rad ** 2
+    # Remove units, get ndarray of values
+    arg = (arg.to(u.dimensionless_unscaled)).value
+
+    return c1 * c2 * exp1(arg)


Test patch:
diff --git a/plasmapy/diagnostics/tests/test_thomson.py b/plasmapy/diagnostics/tests/test_thomson.py
index e16b095e..fbacca44 100644
--- a/plasmapy/diagnostics/tests/test_thomson.py
+++ b/plasmapy/diagnostics/tests/test_thomson.py
@@ -24,10 +24,9 @@ def gen_collective_spectrum():
     """
     wavelengths = np.arange(520, 545, 0.01) * u.nm
     probe_wavelength = 532 * u.nm
-    ne = 5e17 * u.cm ** -3
+    n = 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+"]
@@ -35,10 +34,9 @@ def gen_collective_spectrum():
     alpha, Skw = thomson.spectral_density(
         wavelengths,
         probe_wavelength,
-        ne,
+        n,
         Te,
         Ti,
-        fract=fract,
         ion_species=ion_species,
         probe_vec=probe_vec,
         scatter_vec=scatter_vec,
@@ -55,13 +53,13 @@ def gen_multiple_ion_species_spectrum():
     """
     wavelengths = np.arange(520, 545, 0.01) * u.nm
     probe_wavelength = 532 * u.nm
-    ne = 5e17 * u.cm ** -3
+    n = 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])
+    ifract = 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
+    electron_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
@@ -70,14 +68,14 @@ def gen_multiple_ion_species_spectrum():
     alpha, Skw = thomson.spectral_density(
         wavelengths,
         probe_wavelength,
-        ne,
+        n,
         Te,
         Ti,
-        fract=fract,
+        ifract=ifract,
         ion_species=ion_species,
         probe_vec=probe_vec,
         scatter_vec=scatter_vec,
-        fluid_vel=fluid_vel,
+        electron_vel=electron_vel,
         ion_vel=ion_vel,
     )
 
@@ -91,10 +89,9 @@ def gen_non_collective_spectrum():
     """
     wavelengths = np.arange(500, 570, 0.01) * u.nm
     probe_wavelength = 532 * u.nm
-    ne = 5e15 * u.cm ** -3
+    n = 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+"]
@@ -102,10 +99,9 @@ def gen_non_collective_spectrum():
     alpha, Skw = thomson.spectral_density(
         wavelengths,
         probe_wavelength,
-        ne,
+        n,
         Te,
         Ti,
-        fract=fract,
         ion_species=ion_species,
         probe_vec=probe_vec,
         scatter_vec=scatter_vec,
@@ -119,11 +115,11 @@ 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
+    n = 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
+    ifract = np.array([1.0])
+    Te = np.array([10]) * u.eV
     Ti = np.array([10]) * u.eV
     ion_species = "C-12 5+"
 
@@ -132,10 +128,10 @@ def test_different_input_types():
         alpha, Skw = thomson.spectral_density(
             wavelengths,
             probe_wavelength,
-            ne,
+            n,
             Te,
             Ti,
-            fract=np.array([0.5, 0.5]),
+            ifract=np.array([0.5, 0.5]),
             ion_species=ion_species,
             probe_vec=probe_vec,
             scatter_vec=scatter_vec,
@@ -146,10 +142,10 @@ def test_different_input_types():
         alpha, Skw = thomson.spectral_density(
             wavelengths,
             probe_wavelength,
-            ne,
+            n,
             Te,
             np.array([5, 5]) * u.eV,
-            fract=fract,
+            ifract=ifract,
             ion_species=ion_species,
             probe_vec=probe_vec,
             scatter_vec=scatter_vec,
@@ -160,15 +156,43 @@ def test_different_input_types():
         alpha, Skw = thomson.spectral_density(
             wavelengths,
             probe_wavelength,
-            ne,
+            n,
             Te,
             Ti,
-            fract=fract,
+            ifract=ifract,
             ion_species=[],
             probe_vec=probe_vec,
             scatter_vec=scatter_vec,
         )
 
+    # Raise a Value Error with inconsistent electron array lengths
+    # Te.size != efract.size
+    with pytest.raises(ValueError):
+        alpha, Skw = thomson.spectral_density(
+            wavelengths,
+            probe_wavelength,
+            n,
+            np.array([1, 10]) * u.eV,
+            Ti,
+            efract=np.array([0.5, 0.2, 0.3]),
+            probe_vec=probe_vec,
+            scatter_vec=scatter_vec,
+        )
+
+    # Electron vel shape not compatible with efract.size
+    with pytest.raises(ValueError):
+        alpha, Skw = thomson.spectral_density(
+            wavelengths,
+            probe_wavelength,
+            n,
+            Te,
+            Ti,
+            efract=np.array([0.5, 0.5]),
+            electron_vel=np.array([[100, 0, 0]]) * u.km / u.s,
+            probe_vec=probe_vec,
+            scatter_vec=scatter_vec,
+        )
+
 
 def test_collective_spectrum():
     """
@@ -233,9 +257,9 @@ def test_multiple_ion_species_spectrum():
     max_wavelength = wavelength.value[np.argmax(Skw.value)]
 
     # Check width
-    assert np.isclose(width, 0.1599, 1e-3), (
+    assert np.isclose(width, 0.049999, 1e-2), (
         f"Multiple ion species case spectrum width is {width} instead of "
-        "expected 0.1599"
+        "expected 0.04999"
     )
 
     # Check max value
@@ -245,7 +269,67 @@ def test_multiple_ion_species_spectrum():
     )
 
     # Check max peak location
-    assert np.isclose(max_wavelength, 531.549, 1e-3), (
+    assert np.isclose(max_wavelength, 530.799, 1e-2), (
         "Multiple ion species case spectrum peak wavelength is "
-        f"{max_wavelength} instead of expected 531.549"
+        f"{max_wavelength} instead of expected 530.79"
     )
+
+
+def test_split_populations():
+    """
+    This test makes sure that splitting a single population of ions or electrons
+    into two identical halves returns the same result.
+    """
+
+    wavelengths = np.arange(520, 545, 0.01) * u.nm
+    probe_wavelength = 532 * u.nm
+    n = 5e17 * u.cm ** -3
+    probe_vec = np.array([1, 0, 0])
+    scatter_vec = np.array([0, 1, 0])
+
+    # Combined
+    Te = np.array([10]) * u.eV
+    Ti = np.array([10]) * u.eV
+    ion_species = ["H+"]
+    ifract = np.array([1.0])
+    efract = np.array([1.0])
+
+    alpha, Skw0 = thomson.spectral_density(
+        wavelengths,
+        probe_wavelength,
+        n,
+        Te,
+        Ti,
+        ifract=ifract,
+        efract=efract,
+        ion_species=ion_species,
+        probe_vec=probe_vec,
+        scatter_vec=scatter_vec,
+    )
+
+    # Split e and i populations into two parts
+    # this should not change the results since the parts are identical
+    Te = np.array([10, 10]) * u.eV
+    Ti = np.array([10, 10]) * u.eV
+    ion_species = ["H+", "H+"]
+    ifract = np.array([0.2, 0.8])
+    efract = np.array([0.8, 0.2])
+
+    alpha, Skw1 = thomson.spectral_density(
+        wavelengths,
+        probe_wavelength,
+        n,
+        Te,
+        Ti,
+        ifract=ifract,
+        efract=efract,
+        ion_species=ion_species,
+        probe_vec=probe_vec,
+        scatter_vec=scatter_vec,
+    )
+
+    # Calculate the deviation between the two spectra
+    # (any differences should be in the noise)
+    deviation = (Skw0 - Skw1) / Skw0 * 100
+
+    assert np.all(deviation < 1e-6), "Failed split populations test"
diff --git a/plasmapy/formulary/tests/test_radiation.py b/plasmapy/formulary/tests/test_radiation.py
new file mode 100644
index 00000000..c206ea9f
--- /dev/null
+++ b/plasmapy/formulary/tests/test_radiation.py
@@ -0,0 +1,34 @@
+import astropy.units as u
+import numpy as np
+import pytest
+
+from plasmapy.formulary.radiation import thermal_bremsstrahlung
+from plasmapy.utils.exceptions import PhysicsError
+
+
+def test_thermal_bremsstrahlung():
+    # Test correct spectrum created
+    frequencies = (10 ** np.arange(15, 16, 0.01)) / u.s
+    ne, Te = 1e22 * u.cm ** -3, 1e2 * u.eV
+    ion_species = "H+"
+
+    spectrum = thermal_bremsstrahlung(frequencies, ne, Te, ion_species=ion_species)
+
+    assert np.isclose(np.max(spectrum).value, 128.4, 1), (
+        f"Spectrum maximum is {np.max(spectrum).value} "
+        "instead of expected value 128.4"
+    )
+
+    # Test violates w > wpe limit
+    small_frequencies = (10 ** np.arange(12, 16, 0.01)) / u.s
+    with pytest.raises(PhysicsError):
+        spectrum = thermal_bremsstrahlung(
+            small_frequencies, ne, Te, ion_species=ion_species
+        )
+
+    # Test violates Rayleigh-Jeans limit
+    small_Te = 1 * u.eV
+    with pytest.raises(PhysicsError):
+        spectrum = thermal_bremsstrahlung(
+            frequencies, ne, small_Te, ion_species=ion_species
+        )