CrossSections.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cross sections\n",
    "Plot our elastic and inelastic scattering cross-sections for hydrogen/tritium. Look angular distributions of outgoing particles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import common libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.constants as const"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make font size bigger and allow LaTeX rendering\n",
    "plt.rcParams.update({'font.size': 15})\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rc('font', family='serif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ionisation cross-sections are from M. E. Rudd, *Phys. Rev. A* **44** (1991)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import cross-section models\n",
    "from src.Rudd1991 import RuddXSec\n",
    "from src.Mott import MottXSec\n",
    "import src.Constants as const\n",
    "from src.Stone2002 import ExcitationXSec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Total cross-sections\n",
    "See what our total cross-sections look like across the tritium $\\beta$-decay spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "incidentEArray = np.linspace(1e3, 19e3, 300)  # Incident electron energies in eV\n",
    "mottXSecArray = np.zeros_like(incidentEArray)\n",
    "ionisationXSecArray = np.zeros_like(incidentEArray)\n",
    "for iE, E in enumerate(incidentEArray):\n",
    "    mottXSecArray[iE] = MottXSec(E).TotalXSec()\n",
    "    ionisationXSecArray[iE] = RuddXSec(E).TotalXSec()\n",
    "\n",
    "# Plot the cross-sections\n",
    "plt.figure()\n",
    "plt.plot(incidentEArray / 1e3, mottXSecArray * 1e4, label='Mott')\n",
    "plt.plot(incidentEArray / 1e3, ionisationXSecArray * 1e4, label='Ionisation')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('Incident electron energy [keV]')\n",
    "plt.ylabel('Cross-section [cm$^2$]')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Angular distributions for outgoing electrons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "endpointEnergy = 18.575e3 # eV\n",
    "endpointElasticXSec = MottXSec(endpointEnergy)\n",
    "endpointIonisationXSec = RuddXSec(endpointEnergy)\n",
    "\n",
    "plt.figure(figsize=(14, 5))\n",
    "plt.subplot(1, 2, 1)\n",
    "angles = np.linspace(1e-5, np.pi - 1e-4, 600)\n",
    "plt.plot(angles, endpointElasticXSec.SinglyDifferentialXSec_theta(angles) * 1e4)\n",
    "plt.xlabel('Scattering angle [radians]')\n",
    "plt.yscale('log')\n",
    "plt.ylabel('$d \\sigma / d \\\\rm{cos}\\\\theta$ [cm$^2$]')\n",
    "plt.title(f'Elastic scattering, T = {endpointEnergy / 1e3:.3f} keV')\n",
    "plt.xlim([-0.01, np.pi])\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.plot(angles, endpointIonisationXSec.DoublyDifferentialXSec(endpointEnergy - const.IONIZATIONENERGYH - 1, angles) * 1e4, label='1 eV')\n",
    "plt.plot(angles, endpointIonisationXSec.DoublyDifferentialXSec(endpointEnergy - const.IONIZATIONENERGYH - 5, angles) * 1e4, label='5 eV')\n",
    "plt.plot(angles, endpointIonisationXSec.DoublyDifferentialXSec(endpointEnergy - const.IONIZATIONENERGYH - 20, angles) * 1e4, label='20 eV')\n",
    "plt.plot(angles, endpointIonisationXSec.DoublyDifferentialXSec(endpointEnergy - const.IONIZATIONENERGYH - 50, angles) * 1e4, label='50 eV')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('Scattering angle [radians]')\n",
    "plt.ylabel('$d^2 \\sigma / (d \\\\rm{cos}\\\\theta dW)$ [cm$^2$ eV$^{-1}$]')\n",
    "plt.xlim([0, np.pi / 4])\n",
    "plt.legend(title = 'Secondary energy')\n",
    "plt.title(f'Ionisation, T = {endpointEnergy / 1e3:.3f} keV')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Excitation cross-sections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the excitation cross-section as function of energy\n",
    "exciteXSec = ExcitationXSec(100.0)\n",
    "eArray2p, xSecArr2p = exciteXSec.XSec2pData()\n",
    "eArray3p, xSecArr3p = exciteXSec.XSec3pData()\n",
    "eArray4p, xSecArr4p = exciteXSec.XSec4pData()\n",
    "eArray5p, xSecArr5p = exciteXSec.XSec5pData()\n",
    "eArray6p, xSecArr6p = exciteXSec.XSec6pData()\n",
    "\n",
    "totalXSecArr = np.array([ExcitationXSec(energy).TotalXSec() for energy in eArray2p])\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(eArray2p, xSecArr2p, '-', label='1s-2p')\n",
    "plt.plot(eArray3p, xSecArr3p, '-', label='1s-3p')\n",
    "plt.plot(eArray4p, xSecArr4p, '-', label='1s-4p')\n",
    "plt.plot(eArray5p, xSecArr5p, '-', label='1s-5p')\n",
    "plt.plot(eArray6p, xSecArr6p, '-', label='1s-6p')\n",
    "plt.plot(eArray2p, totalXSecArr, '--', label='Total')\n",
    "plt.xscale('log')\n",
    "plt.ylabel('Cross-section [m$^2$]')\n",
    "plt.xlabel('Energy [eV]')\n",
    "plt.title(r'$\\sigma_{BE}$ for H')\n",
    "plt.legend()\n",
    "plt.ylim([0., 1.1 * np.max(totalXSecArr)])\n",
    "plt.xlim([10., 3e3])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of excitation and ionisation cross-sections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare excitation and ionisation cross-sections\n",
    "energyArray = np.logspace(np.log10(14.0), np.log10(3e3), 300)\n",
    "exciteXSecArray = np.array([ExcitationXSec(E).TotalXSec() for E in energyArray])\n",
    "ionisationXSecArray = np.array([RuddXSec(E).TotalXSec() for E in energyArray])\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(energyArray, exciteXSecArray, label='Excitation')\n",
    "plt.plot(energyArray, ionisationXSecArray, label='Ionisation')\n",
    "plt.xscale('log')\n",
    "plt.ylabel('Cross-section [m$^2$]')\n",
    "plt.xlabel('Energy [eV]')\n",
    "plt.legend()\n",
    "plt.xlim([10., 3e3])\n",
    "plt.ylim([0., 1.1 * np.max(exciteXSecArray)])\n",
    "plt.title('Cross-sections for atomic H')\n",
    "plt.show()"
   ]
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}