{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Emission of Thermal Bremsstrahlung by a Maxwellian Plasma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[thermal_bremsstrahlung]: ../../api/plasmapy.formulary.radiation.thermal_bremsstrahlung.rst#plasmapy.formulary.radiation.thermal_bremsstrahlung\n",
    "\n",
    "The [radiation.thermal_bremsstrahlung][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",
    "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, 0.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
}