NDmodule/ND03_DopplerBroaden.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Doppler broadening of resonances\n",
    "\n",
    "Please indicate your name below, since you will need to submit this notebook.\n",
    "\n",
    "Don't forget to save your progress during the datalab to avoid any loss due to crashes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "name=''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Previously we studied how energy-dependent cross sections can be reconstructed from the resonance parameters. These resonance parameters are a property of the nucleus, hence they do not depend on temperature. Therefore it can be considered that our reconstuction was related to 0 K temperature. In this datalab, we will study how temperature affects the cross sections, and we will perform an approximation to broaden the resonance of U-238 at 6.67 eV. The student is expected to read the chapter \"Effects of nuclear motion\" from [Nuclear Reactor Analysis from Duderstadt and Hamilton](https://deepblue.lib.umich.edu/bitstream/handle/2027.42/89079/1976_Nuclear_Reactor_Analysis.pdf) (page 45-52) along this section. The datalab provides only a brief review of the theory.\n",
    "\n",
    "In the ENDF files the resonance parameters belong to the case when the target nuclei are at rest, however in practice this is not the case due to thermal motion with a thermal energy of $E\\cong kT$. At temperatures typical for nuclear reactors this energy is only comparable to the neutron's kinetic energy if the neutron has also thermal energies. Execute the code cell below to inspect thermal energies at typical reactor temperatures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k=8.617333262145e-5 #Boltzmann constant in eV/K\n",
    "Temps=[293,300,600,900,1200] #in K\n",
    "m=1.67492749804E-27 #kg\n",
    "eVtoJ=1.60217662E-19\n",
    "\n",
    "for T in Temps:\n",
    "    print('Thermal energy of nuclei at {} K: {:.4f} eV; and neutron speed: {:.2f} km/s'.format(T,\n",
    "                                                                               k*T, \n",
    "                                                                               (2*k*T*eVtoJ/m)**(1/2)/1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In later datalabs we will study, how this thermal motion can result in upscattering of low-energy neutrons (ie. neutrons might gain energy in a scattering event).\n",
    "\n",
    "However thermal motion must be taken into account even if the neutron's kinetic energy is way above the thermal energy of the nuclei: namely at resonances, since their width is is comparable with these energies, therefore the energy dependence of the cross sections close to resonances strongly depends on temperatue. This is referred to as Doppler effect or Doppler broadening, and is the subject of this notebook. \n",
    "\n",
    "In case of thermal motion, one needs to take into account the relative speed ($v_r=|\\mathbf{v}-\\mathbf{V}|$) between the neutron and the nuclei. The interaction frequency can be given as\n",
    "\n",
    "$$|\\mathbf{v}-\\mathbf{V}|\\sigma(|\\mathbf{v}-\\mathbf{V}|)N$$\n",
    "\n",
    "where $N$ is the number density of the nuclei.\n",
    "\n",
    "If we assume that the distribution of nuclear velocities is given by the Maxwell-Boltzmann distribution $M(\\mathbf{V},T)$, then one can arrive to the formula describing the thermally averaged cross sections:\n",
    "\n",
    "$$\\bar\\sigma(v,T)=\\frac{1}{v}\\int d^3 V |\\mathbf{v}-\\mathbf{V}|\\sigma(|\\mathbf{v}-\\mathbf{V}|)M(\\mathbf{V},T)$$\n",
    "\n",
    "We can recall from the earlier datalab that the Single-level Breit-Wigner (SLBW) formula for the radiative capture cross section is\n",
    "\n",
    "$$\\sigma_\\gamma(E)=\\sigma_0\\frac{\\Gamma_\\gamma}{\\Gamma}\\Big(\\frac{E_0}{E}\\Big)^{1/2}\\frac{1}{1+y^2}$$\n",
    "\n",
    "with\n",
    "\n",
    "$$y=\\frac{2}{\\Gamma}(E-E_0)$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\\sigma_0=2.608\\cdot 10^6\\frac{(A+1)^2}{A^2E_0(eV)}\\frac{\\Gamma_n}{\\Gamma}g$$\n",
    "\n",
    "\n",
    "So now we \"only\" need to bring everything together to arrive to a formula describing the temperature dependent cross sections. Again, for the proper math we refer the student to the DH book.\n",
    "\n",
    "Let us first reconstruct the first resonance of U-238 with the SLBW formalism as we did during the previous datalab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "E0s, Js, GNs, GGs, GFAs, GFBs = np.loadtxt('data/U238res.dat').transpose() #transpose so we can read the columns\n",
    "E0,J,GN,GG=E0s[0],Js[0],GNs[0],GGs[0]\n",
    "\n",
    "I=0.0\n",
    "A=238\n",
    "G=GN+GG                  #total line width\n",
    "g = (2*J+1)/(2*(2*I+1))  #statistical spin factor\n",
    "sigma_0 = 2.608e6 * (A+1)**2/(A**2 * E0) * (GN/G) * g\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E=np.linspace(6., 7.347, 1000)\n",
    "y = (2/G)*(E - E0)\n",
    "sigma_g_0K = sigma_0 * (GG/G) * np.sqrt(E0/E) * (1/(1+y**2))\n",
    "\n",
    "plt.figure(0)\n",
    "plt.plot(E, sigma_g_0K,'k')\n",
    "\n",
    "plt.xlabel('Energy (eV)')\n",
    "plt.ylabel('Capture cross section (b)')            \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Doppler broadening\n",
    "\n",
    "Now we will try to reconstruct the resonance with the Doppler broadening function. After substituting the SLBW formalism and the Maxwell-Boltzmann distribution into the definition of the thermally averaged cross section, we can arrive to\n",
    "\n",
    "$$\\sigma_\\gamma(E,T)=\\sigma_0\\frac{\\Gamma_\\gamma}{\\Gamma}\\Big(\\frac{E_0}{E}\\Big)^{1/2}\\Psi(\\zeta,x)$$\n",
    "\n",
    "with\n",
    "\n",
    "$$x=2\\frac{E-E0}{\\Gamma}$$\n",
    "\n",
    "$$\\zeta=\\frac{\\Gamma}{\\Gamma_D}$$\n",
    "\n",
    "and the Doppler width is defined as:\n",
    "\n",
    "$$\\Gamma_D=\\Big(\\frac{4E_0kT}{A}\\Big)^{1/2}$$\n",
    "\n",
    "The broadening function is\n",
    "\n",
    "$$\\Psi(\\zeta,x)=\\frac{\\zeta}{2}\\int_{-2E/\\Gamma}^\\infty \\frac{dy}{1+y^2}\\Bigg[\\exp\\Big(-\\frac{(v-v_r)^2}{2v_{th}^2}\\Big)-\\exp\\Big(-\\frac{(v+v_r)^2}{2v_{th}^2}\\Big)\\Bigg]$$\n",
    "\n",
    "This integral can be evaluated numerically. However for simplicity, in the following we will use an approximate expression from Bethe and Placzek (the rationale and the limitation of this approximation is described in detail in the DH book, nevertheless it is rather accurate at lower temperatures and for not too low-energy resonances).\n",
    "\n",
    "$$\\Psi(\\zeta,x) \\sim \\psi(\\zeta,x)=\\frac{\\zeta}{2\\sqrt{\\pi}}\\int_{-\\infty}^{\\infty}dy\\frac{\\exp\\Big[-\\frac{1}{4}(x-y)^2\\zeta^2\\Big]}{1+y^2}$$\n",
    "\n",
    "To make your life easier, this function was implemented as the `psi()` Python function. It uses the `scipy.integrate.quad` method to perform the integration. Excecute and review the code cell below. You can test it, by comparing the returned values for various input values to the tabulated values found in Table 2-2 of the DH book."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import quad\n",
    "\n",
    "def psi(zeta,x):\n",
    "    \"\"\"\n",
    "    Evaluating the Doppler Broadening Psi function (the simplified given by Bethe and Placzek).\n",
    "    See DH Eq. (2.97)\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    zeta : float\n",
    "        zeta parameter, where zeta=Gamma/Gamma_D\n",
    "    x : list or float\n",
    "        x=2*(E-E0)/Gamma. Can be either a single value or a list of values. In case it is a list\n",
    "        psi is evaluated for each x_i value.\n",
    "        \n",
    "    Examples\n",
    "    --------\n",
    "    >>> psi(0.15,4)\n",
    "    0.11268361165162198\n",
    "    >>> psi(0.15,[4,6])\n",
    "    array([0.11268361, 0.10165318])\n",
    "    \n",
    "    Note\n",
    "    ----\n",
    "    The function tests whether it is possible to loop through the values of x. In that case it\n",
    "    returns an array of values. Otherwise it will consider x being an int or float, and returns\n",
    "    a single value. Nevertheless, no further exception is handled (eg the function breaks \n",
    "    if x is a string)\n",
    "    \"\"\"\n",
    "    def integrand(y,zeta,x):\n",
    "        return (zeta/(2*np.sqrt(np.pi)))*(np.exp(-0.25*(x-y)**2*zeta**2))/(1+y**2)\n",
    "    \n",
    "    try:\n",
    "        resu=[]\n",
    "        for xi in x:\n",
    "            resu.append(quad(integrand, -np.inf, np.inf, args=(zeta, xi))[0])\n",
    "        return np.array(resu)\n",
    "    except TypeError:\n",
    "        return quad(integrand, -np.inf, np.inf, args=(zeta, x))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the only thing left is to reconstruct the broadened resonance at various temperatures. Complete the code cell below. And conclude your findings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(E,sigma_g_0K,label='T=0 K')\n",
    "print('The integral under the resonance at T=0 K is {}'.format(np.trapz(sigma_g_0K,E)))\n",
    "\n",
    "for T in [300,600,900,1200]:\n",
    "\n",
    "    GD=#Doppler width \n",
    "    zeta=#zeta parametr\n",
    "    x=#x values\n",
    "    sigma_g =#broadened cross section\n",
    "\n",
    "    plt.plot(E,sigma_g,label='T={} K'.format(T))\n",
    "    print('The integral under the resonance at T={} K is {}'.format(T,np.trapz(sigma_g,E)))\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your conclusion comes here"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}