Repository: PlasmaPy/PlasmaPy
Pull request #892 (2020-08-25 14:30:36)

Function to calculate the thermal bremsstrahlung spectrum
I wrote a function to calculate the power spectrum of thermal bremsstrahlung emitted by a Maxwellian plasma using the equations given on [the Wikipedia page](https://en.wikipedia.org/wiki/Bremsstrahlung#In_plasma). I also created the appropriate documentation (docstring, an ipynb example). 

I wasn't sure where to put this: it doesn't fit into any of the existing Formulary categories so I created a new one for "Radiation". Alternate names could maybe be "EMradiation" or something?

codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/892?src=pr&el=h1) Report
> Merging [#892](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/892?src=pr&el=desc) into [master](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/37120c7cf0c0703b3363e17d678746c6d9ed6f76?el=desc) will **increase** coverage by `0.02%`.
> The diff coverage is `100.00%`.

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

```diff
@@            Coverage Diff             @@
##           master     #892      +/-   ##
==========================================
+ Coverage   96.18%   96.20%   +0.02%     
==========================================
  Files          59       60       +1     
  Lines        5367     5397      +30     
==========================================
+ Hits         5162     5192      +30     
  Misses        205      205              
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/892?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [plasmapy/formulary/radiation.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/892/diff?src=pr&el=tree#diff-cGxhc21hcHkvZm9ybXVsYXJ5L3JhZGlhdGlvbi5weQ==) | `100.00% <100.00%> (ø)` | |

------

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

StanczakDominik: The py36_minimal_conda [linux] failure is unrelated (#894 attempted a fix), so let's merge it :tada:  Thanks!
StanczakDominik: Wrong button!

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/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/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/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
+        )