Jupyter notebook for plotting 3D Compton+Rayleigh scatter intensity

readme.md

A Jupyter notebook for plotting the 3D Compton+Rayleigh scatter intensity for a specified energy and compound

Installation instructions

• Install miniconda or anaconda

• Copy the accompanying files into a folder (call it path)

• In a terminal/shell

  > cd path
  > conda env create -n compton environment.yml
  

Running

  > conda activate compton
  (compton) PS C:\path> jupyter lab

• In Jupyter, open misc_tomopy_mayavi_compton-energy_simplified.ipynb

Note: The mayavi plots will open in their own separate window

environment.yml

name: compton
channels:
  - conda-forge
  - defaults
dependencies:
  - xraylib
  - mayavi
  - jupyterlab
  - matplotlib
  - scipy

misc_tomopy_mayavi_compton-energy_simplified.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculations of Compton scattering and fluorescence into detectors perpendicular to beam axis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Useful resources  \n",
    "The API https://github.com/tschoonj/xraylib/wiki/The-xraylib-API-list-of-all-functions  \n",
    "For some Python usage examples see C:\\Program Files (x86)\\xraylib\\Example\\xrlexample5.py  \n",
    "xraylib The definitive manual  \n",
    "The paper http://dx.doi.org/10.1016/j.sab.2011.09.011"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coordinate system used by DCSP_Rayl and DCSP_Compt methods:\n",
    "\"Given an energy E, a scattering polar angle $\\theta$ and scattering azimuthal angle $\\phi$, returns the Klein-Nishina differential cross section for a polarized beam\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](http://upload.wikimedia.org/wikipedia/commons/thumb/4/4f/3D_Spherical.svg/200px-3D_Spherical.svg.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import xraylib as xrl\n",
    "import mayavi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may need to increase the following value if the want to plot crosssections for lower energies. This seems OK for energies down to ~2keV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "THETA_MIN = 0.015"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize_scalar\n",
    "from mayavi import mlab\n",
    "\n",
    "\n",
    "# From https://stackoverflow.com/questions/36816537\n",
    "def plotsurf(f, show=True, *args, **kwargs):\n",
    "    '''Plot the polar coordinatised function f.\n",
    "\n",
    "    Args:\n",
    "        f: function f(theta, phi) in polar coordinates theta in [0, pi], phi in [0, 2 pi]\n",
    "    \n",
    "    '''\n",
    "    theta, phi = np.linspace(THETA_MIN, np.pi, 80), np.linspace(0, 2 * np.pi, 40)\n",
    "    THETA, PHI = np.meshgrid(theta, phi)\n",
    "    R = f(THETA, PHI)\n",
    "    X = R * np.sin(THETA) * np.cos(PHI)\n",
    "    Y = R * np.sin(THETA) * np.sin(PHI)\n",
    "    Z = R * np.cos(THETA)\n",
    "    \n",
    "    s = mlab.mesh(X, Y, Z)\n",
    "    mlab.show()\n",
    "\n",
    "\n",
    "# Following http://gael-varoquaux.info/programming/mayavi-representing-an-additional-scalar-on-surfaces.html\n",
    "def plotsurf_with_colour(f, value, show=True, *args, **kwargs):\n",
    "    '''Plot the polar coordinatised function f.\n",
    "\n",
    "    Args:\n",
    "        f: radius function f(theta, phi) in polar coordinates theta in [0, pi], phi in [0, 2 pi]\n",
    "        value: scalar function f(theta, phi) to map to colour\n",
    "    \n",
    "    '''\n",
    "    theta, phi = np.linspace(THETA_MIN, np.pi, 150), np.linspace(0, 2 * np.pi, 40)\n",
    "    THETA, PHI = np.meshgrid(theta, phi)\n",
    "    R = f(THETA, PHI)\n",
    "    X = R * np.sin(THETA) * np.cos(PHI)\n",
    "    Y = R * np.sin(THETA) * np.sin(PHI)\n",
    "    Z = R * np.cos(THETA)\n",
    "\n",
    "    scalars = value(THETA, PHI)\n",
    "    assert R.shape == scalars.shape\n",
    "    s = mlab.mesh(X, Y, Z, scalars=scalars)\n",
    "    mlab.colorbar(orientation='vertical', nb_labels=4)\n",
    "\n",
    "    mlab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the incident energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "double ComptonEnergy(double E0, double theta)  \n",
    "double DCSP_Compt(int Z, double E, double theta, double phi)  \n",
    "double DCSP_KN(double E, double theta, double phi)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "vDCSP_KN = np.vectorize(xrl.DCSP_KN, excluded={0})\n",
    "vDCS_KN = np.vectorize(xrl.DCS_KN, excluded={0})\n",
    "vDCSP_Compt_CP = np.vectorize(xrl.DCSP_Compt_CP, excluded={0,1})\n",
    "vDCSP_Rayl_CP = np.vectorize(xrl.DCSP_Rayl_CP, excluded={0,1})\n",
    "vComptonEnergy = np.vectorize(xrl.ComptonEnergy, excluded={0,2})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "E0 = 100.0    # Energy (keV)\n",
    "f = lambda t, p: vDCSP_KN(E0, t, p)\n",
    "value = lambda t, p: vComptonEnergy(E0, t)     # Surface colouring function\n",
    "plotsurf_with_colour(f, value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "E0 = 30.0\n",
    "f = lambda t, p: vDCS_KN(E0, t)\n",
    "plotsurf(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "    DCSP_Compt_CP(char const [] compound, double E, double theta, double phi) -> double\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    compound: char const []\n",
      "    E: double\n",
      "    theta: double\n",
      "    phi: double\n",
      "\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "print(xrl.DCSP_Compt_CP.__doc__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NIST Compounds understood by xraylib are https://github.com/tschoonj/xraylib/blob/master/src/xraylib-nist-compounds-internal.h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "E0 = 30.0\n",
    "material = 'Al'\n",
    "# material = 'Polymethyl Methacralate (Lucite, Perspex)'\n",
    "f = lambda t, p: vDCSP_Rayl_CP(material, E0, t, p)\n",
    "plotsurf(f, alpha=0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "E0 = 3.0\n",
    "# material = 'Al'\n",
    "# material = 'Polymethyl Methacralate (Lucite, Perspex)'\n",
    "material = 'Air, Dry (near sea level)'\n",
    "f = lambda t, p: vDCSP_Compt_CP(material, E0, t, p)\n",
    "energy = lambda t, p: vComptonEnergy(E0, t)\n",
    "plotsurf_with_colour(f, energy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "E0 = 2.0\n",
    "# material = 'Al'\n",
    "material = 'Polymethyl Methacralate (Lucite, Perspex)'\n",
    "# material = 'Air, Dry (near sea level)'\n",
    "f = lambda t, p: vDCSP_Compt_CP(material, E0, t, p) + vDCSP_Rayl_CP(material, E0, t, p)\n",
    "energy = lambda t, p: vComptonEnergy(E0, t)\n",
    "plotsurf_with_colour(f, energy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:vis]",
   "language": "python",
   "name": "conda-env-vis-py"
  },
  "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
}