Merge pull request #497 from scipp/scattering-table

Add a function to access scattering parameters

docs/api-reference/index.md

@@ -68,6 +68,7 @@ and possible confusion of `theta` (from Bagg's law) with `theta` in spherical co
    :template: module-template.rst
    :recursive:
 
+   atoms
    conversion
    io
    logging

docs/bibliography.bib

@@ -20,3 +20,15 @@ @article{busing:1967
     volume = {22},
     pages = {457--464},
 }
+
+@article{sears:1992,
+    author = {Varley F. Sears},
+    title = {Neutron scattering lengths and cross sections},
+    journal = {Neutron News},
+    volume = {3},
+    number = {3},
+    pages = {26-37},
+    year = {1992},
+    publisher = {Taylor & Francis},
+    doi = {10.1080/10448639208218770},
+}

pyproject.toml

@@ -95,7 +95,7 @@ exclude_dirs = ["docs/conf.py", "tests", "tools"]
 
 [tool.codespell]
 ignore-words-list = "elemt"
-skip = "./.git,./.tox,*/.virtual_documents,*/.ipynb_checkpoints,*.pdf,*.svg"
+skip = "./.git,./.tox,*/.virtual_documents,*/.ipynb_checkpoints,*.pdf,*.svg,*.csv"
 
 [tool.black]
 skip-string-normalization = true

src/scippneutron/__init__.py

@@ -41,6 +41,7 @@
 from .instrument_view import instrument_view
 from .io.nexus.load_nexus import load_nexus, load_nexus_json
 from .data_streaming.data_stream import data_stream
+from . import atoms
 from . import data
 
 del importlib

src/scippneutron/atoms/__init__.py

@@ -0,0 +1,130 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
+"""Parameters for neutron interactions with atoms."""
+from __future__ import annotations
+
+import dataclasses
+import importlib.resources
+from functools import lru_cache
+from typing import Optional, TextIO, Union
+
+import scipp as sc
+
+
+def reference_wavelength() -> sc.Variable:
+    """Return the reference wavelength for absorption cross-sections.
+
+    Returns
+    -------
+    :
+        1.7982 Å
+    """
+    return sc.scalar(1.7982, unit='angstrom')
+
+
+@dataclasses.dataclass(frozen=True, eq=False)
+class ScatteringParams:
+    """Scattering parameters for neutrons with a specific element / isotope.
+
+    Provides access to the scattering lengths and cross-sections of neutrons
+    with a given element or isotope.
+    Values have been retrieved at 2024-02-19T17:00:00Z from the list at
+    https://www.ncnr.nist.gov/resources/n-lengths/list.html
+    which is based on :cite:`sears:1992`.
+    Values are ``None`` where the table does not provide values.
+
+    The absorption cross-section applies to neutrons with a wavelength
+    of 1.7982 Å.
+    See :func:`reference_wavelength`.
+    """
+
+    isotope: str
+    """Element / isotope name."""
+    coherent_scattering_length_re: Optional[sc.Variable]
+    """Bound coherent scattering length (real part)."""
+    coherent_scattering_length_im: Optional[sc.Variable]
+    """Bound coherent scattering length (imaginary part)."""
+    incoherent_scattering_length_re: Optional[sc.Variable]
+    """Bound incoherent scattering length (real part)."""
+    incoherent_scattering_length_im: Optional[sc.Variable]
+    """Bound incoherent scattering length (imaginary part)."""
+    coherent_scattering_cross_section: Optional[sc.Variable]
+    """Bound coherent scattering cross-section."""
+    incoherent_scattering_cross_section: Optional[sc.Variable]
+    """Bound incoherent scattering cross-section."""
+    total_scattering_cross_section: Optional[sc.Variable]
+    """Total bound scattering cross-section."""
+    absorption_cross_section: Optional[sc.Variable]
+    """Absorption cross-section for λ = 1.7982 Å neutrons."""
+
+    def __eq__(self, other: object) -> Union[bool, type(NotImplemented)]:
+        if not isinstance(other, ScatteringParams):
+            return NotImplemented
+        return all(
+            self.isotope == other.isotope
+            if field.name == 'isotope'
+            else _eq_or_identical(getattr(self, field.name), getattr(other, field.name))
+            for field in dataclasses.fields(self)
+        )
+
+    @staticmethod
+    @lru_cache()
+    def for_isotope(isotope: str) -> ScatteringParams:
+        """Return the scattering parameters for the given element / isotope.
+
+        Parameters
+        ----------
+        isotope:
+            Name of the element or isotope.
+            For example, 'H', '3He', 'V', '50V'.
+
+        Returns
+        -------
+        :
+            Neutron scattering parameters.
+        """
+        with _open_scattering_parameters_file() as f:
+            while line := f.readline():
+                name, rest = line.split(',', 1)
+                if name == isotope:
+                    return _parse_line(isotope, rest)
+        raise ValueError(f"No entry for element / isotope '{isotope}'")
+
+
+def _open_scattering_parameters_file() -> TextIO:
+    return (
+        importlib.resources.files('scippneutron.atoms')
+        .joinpath('scattering_parameters.csv')
+        .open('r')
+    )
+
+
+def _parse_line(isotope: str, line: str) -> ScatteringParams:
+    line = line.rstrip().split(',')
+    return ScatteringParams(
+        isotope=isotope,
+        coherent_scattering_length_re=_assemble_scalar(line[0], line[1], 'fm'),
+        coherent_scattering_length_im=_assemble_scalar(line[2], line[3], 'fm'),
+        incoherent_scattering_length_re=_assemble_scalar(line[4], line[5], 'fm'),
+        incoherent_scattering_length_im=_assemble_scalar(line[6], line[7], 'fm'),
+        coherent_scattering_cross_section=_assemble_scalar(line[8], line[9], 'barn'),
+        incoherent_scattering_cross_section=_assemble_scalar(
+            line[10], line[11], 'barn'
+        ),
+        total_scattering_cross_section=_assemble_scalar(line[12], line[13], 'barn'),
+        absorption_cross_section=_assemble_scalar(line[14], line[15], 'barn'),
+    )
+
+
+def _assemble_scalar(value: str, std: str, unit: str) -> Optional[sc.Variable]:
+    if not value:
+        return None
+    value = float(value)
+    variance = float(std) ** 2 if std else None
+    return sc.scalar(value, variance=variance, unit=unit)
+
+
+def _eq_or_identical(a: Optional[sc.Variable], b: Optional[sc.Variable]) -> bool:
+    if a is None:
+        return b is None
+    return sc.identical(a, b)