feat: allow computing transmission map at any point (not only infinit…

…ely far away)

docs/user-guide/absorption-correction.ipynb

@@ -0,0 +1,245 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "a6de8431-db39-49a5-b5dc-1afee82bc7d2",
+   "metadata": {},
+   "source": [
+    "# Absorption correction for cylindrical samples\n",
+    "\n",
+    "The intensity at a detector element located at $\\mathbf{x}$ is modeled as\n",
+    "$$\n",
+    "I(\\lambda, \\mathbf{x}) = \\int_{sample} S(\\lambda, \\arccos\\big(\\frac{\\hat{\\mathbf{b}} \\cdot (\\mathbf{x}-\\mathbf{x}')}{||\\mathbf{x}-\\mathbf{x}'||}\\big)) T(\\lambda, \\mathbf{x}, \\mathbf{x}') \\ d\\mathbf{x}'\n",
+    "$$\n",
+    "where $S(\\lambda, \\theta)$ is the probability that a neutron with wavelength $\\lambda$ scatters with the scattering angle $\\theta$ and $T(\\lambda, \\mathbf{x}, \\mathbf{x}')$ is the transmission rate for neutrons scattering at $\\mathbf{x}'$ towards $\\mathbf{x}$, and $\\hat{\\mathbf{b}}$ is the incoming beam direction.\n",
+    "\n",
+    "If the detector is far away from the sample relative to the size of the sample, then $\\frac{\\hat{\\mathbf{b}} \\cdot (\\mathbf{x}-\\mathbf{x}')}{||\\mathbf{x}-\\mathbf{x}'||}$ is close to constant and the $S$ term can be moved outside of the integtral:\n",
+    "$$\n",
+    "I(\\lambda, \\mathbf{x}) \\approx S(\\lambda, \\theta(\\mathbf{x})) \\int_{sample}  T(\\lambda, \\mathbf{x}, \\mathbf{x}') \\ d\\mathbf{x}'.\n",
+    "$$\n",
+    "To compute the scattering probabiltiy $S$ from the intensity $I$ we need to estimate the term\n",
+    "$$\n",
+    "C(\\lambda, \\mathbf{x}) = \\int_{sample}  T(\\lambda, \\mathbf{x}, \\mathbf{x}') \\ d\\mathbf{x}'.\n",
+    "$$\n",
+    "This is the \"absorption correction\" term.\n",
+    "\n",
+    "The transmission fraction is a function of path length $L$ of the neutron going through the sample\n",
+    "$$\n",
+    "T(\\lambda, \\mathbf{x}, \\mathbf{x}') = \\exp{(-\\mu(\\lambda) L(\\mathbf{x}, \\mathbf{x}'))}\n",
+    "$$\n",
+    "where $\\mu$ is material dependent.\n",
+    "\n",
+    "The path length through the sample depends on the shape of the sample, the scattering point $\\mathbf{x}'$ and the detection position $\\mathbf{x}$.\n",
+    "\n",
+    "To compute $C(\\lambda, \\mathbf{x})$ you can use\n",
+    "```python\n",
+    "scippneutron.absorption.base.compute_transmission_map(\n",
+    "    shape,\n",
+    "    material,\n",
+    "    beam_direction,\n",
+    "    wavelength,\n",
+    "    detector_position\n",
+    ")\n",
+    "```\n",
+    "where `shape` and `material` are sample properties, `beam_direction` corresponds to $\\mathbf{b}$, `wavelength` corresponds to $\\lambda$ and `detector_position` corresponds to $\\mathbf{x}$. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "07bd8cdb-c18a-48e4-a60a-041a0d4cd07c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import scipp as sc\n",
+    "import numpy as np\n",
+    "\n",
+    "from scippneutron.absorption import compute_transmission_map\n",
+    "from scippneutron.absorption.cylinder import Cylinder\n",
+    "\n",
+    "\n",
+    "class DummyMaterial:\n",
+    "    'This is probably not a good physical model, but it illustrates the interface'\n",
+    "    def absorption_coefficient(wavelength):\n",
+    "        return sc.scalar(1, unit='1/(cm * angstrom)') * wavelength\n",
+    "\n",
+    "\n",
+    "sample_shape = Cylinder(\n",
+    "    symmetry_line=sc.vector([0, 1, 0]),\n",
+    "    center_of_base=sc.vector([0, -2, 0], unit='cm'),\n",
+    "    radius=sc.scalar(1, unit='cm'),\n",
+    "    height=sc.scalar(4, unit='cm')\n",
+    ")\n",
+    "\n",
+    "# Create some dummy detector positions.\n",
+    "\n",
+    "## Note\n",
+    "## If the detector array has a large number of elements\n",
+    "## it's best for performance to not evaluate the correction for every detector element\n",
+    "## but instead to evaluate it over some grid of positions and\n",
+    "## interpolate to the positions of the detector elements.\n",
+    "theta = sc.linspace('theta', 0, np.pi, 100, endpoint=True, unit='rad')\n",
+    "phi = sc.linspace('phi', 0, 2 * np.pi, 200, endpoint=False, unit='rad')\n",
+    "r = sc.array(dims='r', values=[2.], unit='m')\n",
+    "dims = (*r.dims, *theta.dims, *phi.dims)\n",
+    "\n",
+    "# Detector positions in grid on a sphere around the origin\n",
+    "detector_positions = sc.vectors(\n",
+    "    dims=dims,\n",
+    "    values=sc.concat(\n",
+    "        [\n",
+    "            r * sc.sin(theta) * sc.cos(phi),\n",
+    "            r * sc.sin(theta) * sc.sin(phi),\n",
+    "            sc.broadcast(r * sc.cos(theta), sizes={**r.sizes, **theta.sizes, **phi.sizes}),\n",
+    "        ],\n",
+    "        dim='row',\n",
+    "    )\n",
+    "    .transpose([*dims, 'row'])\n",
+    "    .values,\n",
+    "    unit=r.unit,\n",
+    ")\n",
+    "\n",
+    "def transmission(quadrature_kind):\n",
+    "    'Evaluates the transmission map for different kinds of quadrature'\n",
+    "    da = compute_transmission_map(\n",
+    "        sample_shape,\n",
+    "        DummyMaterial,\n",
+    "        # Assume beam goes in z-direction\n",
+    "        beam_direction=sc.vector([0, 0, 1]),\n",
+    "        wavelength=sc.linspace('wavelength', 0.5, 2.5, 10, unit='angstrom'),\n",
+    "        detector_position=detector_positions,\n",
+    "        quadrature_kind=quadrature_kind,\n",
+    "    )\n",
+    "    da.coords['phi'] = phi\n",
+    "    da.coords['theta'] = theta\n",
+    "    return da\n",
+    "\n",
+    "\n",
+    "def show_correction_map(da):\n",
+    "    'Plotting utility'\n",
+    "    return (\n",
+    "        da['wavelength', 0]['r', 0].plot() /\n",
+    "        da['wavelength', 0]['r', 0]['theta', da.sizes['theta']//2].plot() /\n",
+    "        da['wavelength', 0]['r', 0]['theta', da.sizes['theta']//2]['phi', da.sizes['phi']//4 - da.sizes['phi']//6:da.sizes['phi']//4 + da.sizes['phi']//6].plot() /\n",
+    "        da['wavelength', 0]['r', 0]['theta', da.sizes['theta']//2]['phi', da.sizes['phi']//2 - da.sizes['phi']//6:da.sizes['phi']//2 + da.sizes['phi']//6].plot() /\n",
+    "        da['wavelength', 0]['r', 0]['phi', da.sizes['phi']//2].plot()\n",
+    "    )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6f03111b-8d1c-4c82-a9e9-b2fcde5c1003",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "transmission('cheap')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "abb0d783-9a57-4daa-b4dc-7295044beb78",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "show_correction_map(transmission('cheap'))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1512c244-118e-40a0-ac88-1510504603da",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "show_correction_map(transmission('medium'))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d34de43c-0341-4656-b739-656cb3aff4a1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "show_correction_map(transmission('expensive'))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "93a0af6e-dd16-4ee4-838d-df2ef014ff4a",
+   "metadata": {},
+   "source": [
+    "## Error estimate"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dbbdbf77-07b0-491e-b898-6daccc03ccd8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import plopp as pp\n",
+    "cheap = transmission('cheap')['wavelength', 0]['r', 0]\n",
+    "medium = transmission('medium')['wavelength', 0]['r', 0]\n",
+    "expensive = transmission('expensive')['wavelength', 0]['r', 0]\n",
+    "\n",
+    "t = cheap.sizes['theta'] // 2\n",
+    "pmin = cheap.sizes['phi'] // 2 - cheap.sizes['phi'] // 8\n",
+    "pmax = cheap.sizes['phi'] // 2 + cheap.sizes['phi'] // 8\n",
+    "pp.plot(\n",
+    "    {\n",
+    "        'cheap': cheap['theta', t]['phi', pmin:pmax],\n",
+    "        'medium': medium['theta', t]['phi', pmin:pmax],\n",
+    "        'expensive': expensive['theta', t]['phi', pmin:pmax],\n",
+    "    }\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "180ab825-0df5-4422-91cf-474523c63fe6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Relative difference between cheap and expensive quadrature\n",
+    "(sc.abs(cheap - expensive)/sc.abs(expensive)).mean().value"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5a3fb34e-449b-4a4c-8ea6-849b0c5166cc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Relative difference between medium and expensive quadrature\n",
+    "(sc.abs(medium - expensive)/sc.abs(expensive)).mean().value"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

src/scippneutron/absorption/__init__.py

@@ -0,0 +1,4 @@
+from .base import compute_transmission_map
+
+
+__all__ = (compute_transmission_map,)