Skip to content

Instantly share code, notes, and snippets.

@gazzar

gazzar/environment.yml

Created December 19, 2022 04:16
Show Gist options
  • Save gazzar/d625e39285829f15ab17337b1411daf2 to your computer and use it in GitHub Desktop.
Save gazzar/d625e39285829f15ab17337b1411daf2 to your computer and use it in GitHub Desktop.
Download ZIP
Jupyter notebook for plotting 3D Compton+Rayleigh scatter intensity
Raw
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

Raw
environment.yml
name: compton
channels:
- conda-forge
- defaults
dependencies:
- xraylib
- mayavi
- jupyterlab
- matplotlib
- scipy
Display the source blob
Display the rendered blob
Raw
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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment