diff --git a/docs/_static/beamline/beamline_coordinates_dark.svg b/docs/_static/beamline/beamline_coordinates_dark.svg
index 905031d64..726fe7ee8 100644
--- a/docs/_static/beamline/beamline_coordinates_dark.svg
+++ b/docs/_static/beamline/beamline_coordinates_dark.svg
@@ -3,438 +3,460 @@
+ id="path147" />
diff --git a/docs/_static/beamline/beamline_coordinates_light.svg b/docs/_static/beamline/beamline_coordinates_light.svg
index 34a9b072c..12200659f 100644
--- a/docs/_static/beamline/beamline_coordinates_light.svg
+++ b/docs/_static/beamline/beamline_coordinates_light.svg
@@ -8,8 +8,29 @@
version="1.1"
id="svg5"
xml:space="preserve"
+ sodipodi:docname="beamline_coordinates_light.svg"
+ inkscape:version="1.3.2 (091e20ef0f, 2023-11-25, custom)"
+ xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
+ xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
xmlns="http://www.w3.org/2000/svg"
- xmlns:svg="http://www.w3.org/2000/svg">gravitygravitypositionposition(detector)sample_positionsource_position(detector)sample_positionsource_positionzzxxyygravityincident_beamgravityincident_beamscattered_beamscattered_beamp
+ y="42.476509">p
diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py
index 65019e3e0..bb624f1e0 100644
--- a/src/scippneutron/conversion/beamline.py
+++ b/src/scippneutron/conversion/beamline.py
@@ -41,16 +41,16 @@
consistently, ``incident_beam``, ``scattered_beam``, and ``two_theta`` can be
computed without knowledge of the coordinate system.
- However, when we need ``phi``, or need to correct ``two_theta`` for gravity,
- we need the actual coordinate system.
+ However, we need the actual coordinate system to compute ``phi``.
ScippNeutron uses a coordinate system aligned with
the incident beam and gravity.
ScippNeutron's coordinate system corresponds to that of
-`NeXus `_.
+`NeXus `_
+when the incident beam is perpendicular to gravity.
The image below shows how coordinates are defined with respect to the
quantities defined above.
-The plot on the right-hand side shows the view from the sample towards the source,
+The plot on the right-hand side shows the view from the sample along the :math:`z`-axis,
that is, the :math:`z`-axis points towards the viewer.
(The sample is placed in the origin here; this is only done for illustration purposes
and not required.)
@@ -71,16 +71,16 @@
.. math::
- \hat{e}_z &= b_1 / |b_1| \\
\hat{e}_y &= -g / |g| \\
+ z_{\text{proj}} &= b_1 - (b_1 \cdot \hat{e}_y) \hat{e}_y \\
+ \hat{e}_z &= z_{\text{proj}} / |z_{\text{proj}}| \\
\hat{e}_x &= \hat{e}_y \times \hat{e}_z
which span an orthogonal, right-handed coordinate system.
Here, :math:`b_1` is the ``incident_beam`` and :math:`g` is the gravity vector.
-This means that the z-axis is parallel to the incident beam and the y-axis is
-antiparallel to gravity.
-Gravity must be orthogonal to the incident beam for this definition to produce
-and orthogonal coordinate system.
+This means that the y-axis is antiparallel to gravity.
+The z-axis is defined by projecting the incident beam onto a plane perpendicular
+to gravity.
Basis vectors can be computed using :func:`beam_aligned_unit_vectors`.
:math:`p = \sqrt{x^2 + y^2}` is the projection of the
@@ -326,7 +326,13 @@ def two_theta(
# operations'.
b1 = incident_beam / L1(incident_beam=incident_beam)
b2 = scattered_beam / L2(scattered_beam=scattered_beam)
- return 2 * sc.atan2(y=sc.norm(b1 - b2), x=sc.norm(b1 + b2))
+
+ y = sc.norm(b1 - b2)
+ b2 += b1
+ x = sc.norm(b2)
+ res = sc.atan2(y=y, x=x, out=x)
+ res *= 2
+ return res
class BeamAlignedUnitVectors(TypedDict):
@@ -349,8 +355,9 @@ def beam_aligned_unit_vectors(
.. math::
- \hat{e}_z &= b_1 / |b_1| \\
\hat{e}_y &= -g / |g| \\
+ z_{\text{proj}} &= b_1 - (b_1 \cdot \hat{e}_y) \hat{e}_y \\
+ \hat{e}_z &= z_{\text{proj}} / |z_{\text{proj}}| \\
\hat{e}_x &= \hat{e}_y \times \hat{e}_z
where :math:`b_1` is the ``incident_beam`` and :math:`g` is the gravity vector.
@@ -367,23 +374,27 @@ def beam_aligned_unit_vectors(
-------
A dict containing the unit vectors with keys 'ex', 'ey', and 'ez'.
+ Raises
+ ------
+ ValueError
+ If the incident beam is parallel to gravity.
+ The rotation of the x-z plane is ill-define din this case.
+
See Also
--------
straight_incident_beam:
Compute the incident beam for a straight beamline.
"""
- if sc.any(
- abs(sc.dot(gravity, incident_beam))
- > sc.scalar(1e-10, unit=incident_beam.unit) * sc.norm(gravity)
- ):
+ ey = -gravity / sc.norm(gravity)
+ # Project incident_beam onto a plane perpendicular to ey.
+ z = incident_beam - sc.dot(incident_beam, ey) * ey
+ z_norm = sc.norm(z)
+ if sc.any(z_norm < sc.scalar(1e-10, unit=z_norm.unit)):
raise ValueError(
- '`gravity` and `incident_beam` must be orthogonal. '
- f'Got a deviation of {sc.dot(gravity, incident_beam).max():c}. '
- 'This is required to fully define spherical coordinates theta and phi.'
+ "Cannot construct a coordinate system. The incident beam and "
+ "gravity are parallel to each other."
)
-
- ez = incident_beam / sc.norm(incident_beam)
- ey = -gravity / sc.norm(gravity)
+ ez = z / sc.norm(z)
ex = sc.cross(ey, ez)
return {
'beam_aligned_unit_x': ex,
@@ -468,18 +479,33 @@ def scattering_angles_with_gravity(
.. math::
x'_d &= x_d \\
- y'_d &= y_d + \frac{|g| m_n^2}{2 h^2} L_2^{\prime\, 2} \lambda^2 \\
+ y'_d &= y_d + \delta_y \\
z'_d &= z_d
Where :math:`|g|` is the strength of gravity, :math:`m_n` is the neutron mass,
- :math:`h` is the Planck constant, and :math:`\lambda` is the wavelength.
+ :math:`h` is the Planck constant, :math:`\lambda` is the wavelength, and
+
+ .. math::
+
+ \delta_y = \frac{|g| m_n^2}{2 h^2} L_2^{\prime\, 2} \lambda^2
+
This gives the gravity-corrected scattering angles:
.. math::
- \mathsf{tan}(2\theta) &= \frac{\sqrt{x_d^2 + y_d^{\prime\, 2}}}{z_d} \\
+ 2\theta &= \sphericalangle(b_1, b_2 + \delta_y \hat{e}_y) \\
\mathsf{tan}(\phi) &= \frac{y'_d}{x_d}
+ where :math:`\sphericalangle` is the angle between two vectors as implemented by
+ :func:`two_theta`.
+ When :math:`b_1` is orthogonal to gravity, we can equivalently use
+
+ .. math::
+
+ \mathsf{tan}(2\theta) = \frac{\sqrt{x_d^2 + y_d^{\prime\, 2}}}{z_d}
+
+ This equation allowed for a more efficient implementation.
+
Attention
---------
The above equation for :math:`y'_d` contains :math:`L_2^{\prime\, 2} = |b'_2|`
@@ -517,6 +543,63 @@ def scattering_angles_with_gravity(
Ignores the ``x`` component when computing ``theta``.
This is used in reflectometry.
"""
+ if sc.any(
+ abs(sc.dot(gravity, incident_beam))
+ > sc.scalar(1e-10, unit=incident_beam.unit) * sc.norm(gravity)
+ ):
+ return _scattering_angles_with_gravity_generic(
+ incident_beam=incident_beam,
+ scattered_beam=scattered_beam,
+ wavelength=wavelength,
+ gravity=gravity,
+ )
+ return _scattering_angles_with_gravity_orthogonal_coords(
+ incident_beam=incident_beam,
+ scattered_beam=scattered_beam,
+ wavelength=wavelength,
+ gravity=gravity,
+ )
+
+
+def _scattering_angles_with_gravity_generic(
+ incident_beam: sc.Variable,
+ scattered_beam: sc.Variable,
+ wavelength: sc.Variable,
+ gravity: sc.Variable,
+) -> SphericalCoordinates:
+ unit_vectors = beam_aligned_unit_vectors(
+ incident_beam=incident_beam, gravity=gravity
+ )
+ ex = unit_vectors['beam_aligned_unit_x']
+ ey = unit_vectors['beam_aligned_unit_y']
+
+ drop_distance = _drop_due_to_gravity(
+ distance=sc.norm(scattered_beam), wavelength=wavelength, gravity=gravity
+ )
+ y = drop_distance + sc.dot(scattered_beam, ey).to(
+ dtype=elem_dtype(wavelength), copy=False
+ )
+ x = sc.dot(scattered_beam, ex).to(dtype=elem_dtype(y), copy=False)
+ phi = sc.atan2(y=y, x=x, out=y)
+
+ drop = drop_distance * (gravity / sc.norm(gravity))
+ drop += scattered_beam
+ return {
+ 'two_theta': two_theta(incident_beam=incident_beam, scattered_beam=drop).to(
+ dtype=elem_dtype(wavelength), copy=False
+ ),
+ 'phi': phi,
+ }
+
+
+# This is an optimized implementation that only works when `incident_beam` and `gravity`
+# are orthogonal to each other. It uses less memory than the generic version.
+def _scattering_angles_with_gravity_orthogonal_coords(
+ incident_beam: sc.Variable,
+ scattered_beam: sc.Variable,
+ wavelength: sc.Variable,
+ gravity: sc.Variable,
+) -> SphericalCoordinates:
unit_vectors = beam_aligned_unit_vectors(
incident_beam=incident_beam, gravity=gravity
)
@@ -609,6 +692,16 @@ def scattering_angle_in_yz_plane(
Includes the ``x`` component when computing ``theta``.
This is used in techniques other than reflectometry.
"""
+ if sc.any(
+ abs(sc.dot(gravity, incident_beam))
+ > sc.scalar(1e-10, unit=incident_beam.unit) * sc.norm(gravity)
+ ):
+ raise ValueError(
+ '`gravity` and `incident_beam` must be orthogonal. '
+ f'Got a deviation of {sc.dot(gravity, incident_beam).max():c}. '
+ 'It is unclear how the angle is defined in this case.'
+ )
+
unit_vectors = beam_aligned_unit_vectors(
incident_beam=incident_beam, gravity=gravity
)
diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py
index 0039974ae..f3d6ac087 100644
--- a/tests/conversion/beamline_conversions_test.py
+++ b/tests/conversion/beamline_conversions_test.py
@@ -7,10 +7,24 @@
import scipp as sc
import scipp.constants
import scipp.testing
+from hypothesis import given, settings
+from hypothesis import strategies as st
+from hypothesis.extra import numpy as npst
from scippneutron.conversion import beamline
+def rotations() -> st.SearchStrategy[sc.Variable]:
+ coefficients = st.floats(
+ allow_nan=False, allow_infinity=False, allow_subnormal=False
+ )
+ return (
+ npst.arrays(dtype=np.float64, shape=4, elements=coefficients)
+ .filter(lambda q: np.linalg.norm(q) != 0)
+ .map(lambda q: sc.spatial.rotation(value=q / np.linalg.norm(q)))
+ )
+
+
def test_straight_incident_beam():
source_position = sc.vectors(
dims=['siti'], values=[[40, 80, 20], [30, 10, 50]], unit='mm'
@@ -157,27 +171,56 @@ def test_two_theta_invariant_under_reflection_about_incident_beam():
sc.testing.assert_allclose(two_theta[0], two_theta[1])
-def test_scattering_angles_requires_gravity_orthogonal_to_incident_beam():
+@given(rotation=rotations())
+@settings(max_examples=10)
+def test_two_theta_invariant_under_rotation(rotation: sc.Variable) -> None:
incident_beam = sc.vector([0.564, 1.2, -10.4], unit='m')
scattered_beam = sc.vectors(
dims=['beam'], values=[[13, 24, 35], [51, -42, 33]], unit='m'
)
+ original = beamline.two_theta(
+ incident_beam=incident_beam, scattered_beam=scattered_beam
+ )
+
+ rotated = beamline.two_theta(
+ incident_beam=rotation * incident_beam, scattered_beam=rotation * scattered_beam
+ )
+
+ sc.testing.assert_allclose(original, rotated)
+
+
+@given(incident_rotation=rotations(), scattered_rotation=rotations())
+@settings(max_examples=50)
+def test_scattering_angles_with_gravity_small_gravity(
+ incident_rotation: sc.Variable, scattered_rotation: sc.Variable
+) -> None:
+ # This case is unphysical but tests the consistency with `two_theta`.
+ incident_beam = sc.vector([0.564, 0.0, 10.4], unit='m')
+ scattered_beam = sc.vectors(
+ dims=['beam'],
+ values=[[13, 24, 35], [51, -42, 33], [4, 23, -17], [-19, -31, 5]],
+ unit='m',
+ )
+ incident_beam = incident_rotation * incident_beam
+ scattered_beam = scattered_rotation * scattered_beam
wavelength = sc.array(dims=['wavelength'], values=[1.2, 1.6, 1.8], unit='Å')
- gravity = sc.vector([0, 0, sc.constants.g.value], unit=sc.constants.g.unit)
+ gravity = sc.vector([0, -1e-11, 0], unit=sc.constants.g.unit)
- with pytest.raises(
- ValueError, match='`gravity` and `incident_beam` must be orthogonal'
- ):
- beamline.scattering_angles_with_gravity(
- incident_beam=incident_beam,
- scattered_beam=scattered_beam,
- wavelength=wavelength,
- gravity=gravity,
- )
+ res = beamline.scattering_angles_with_gravity(
+ incident_beam=incident_beam,
+ scattered_beam=scattered_beam,
+ wavelength=wavelength,
+ gravity=gravity,
+ )
+ expected = beamline.two_theta(
+ incident_beam=incident_beam, scattered_beam=scattered_beam
+ ).broadcast(dims=['wavelength', 'beam'], shape=[3, 4])
+ sc.testing.assert_allclose(res['two_theta'], expected)
-def test_scattering_angles_with_gravity_small_gravity():
+def test_scattering_angles_with_gravity_small_gravity_orthogonal_coords() -> None:
# This case is unphysical but tests the consistency with `two_theta`.
+ # Here, incident_beam and gravity are orthogonal to use the optimized code.
incident_beam = sc.vector([0.564, 0.0, 10.4], unit='m')
scattered_beam = sc.vectors(
dims=['beam'],
@@ -207,6 +250,43 @@ def test_scattering_angles_with_gravity_reproduces_angles(
# This case is unphysical but tests that the function reproduces
# the expected angles using a rotated vector.
+ gravity = sc.vector([1e-12, -1e-11, 2e-12], unit='cm/s^2')
+ incident_beam = sc.vector([0.0, 0.0, 968.0], unit='cm')
+
+ # With this definition, the x-axis has azimuthal=0.
+ rot1 = sc.spatial.rotations_from_rotvecs(sc.vector([-polar, 0, 0], unit='rad'))
+ rot2 = sc.spatial.rotations_from_rotvecs(
+ sc.vector([0, 0, azimuthal - np.pi / 2], unit='rad')
+ )
+ scattered_beam = rot2 * (rot1 * incident_beam)
+
+ wavelength = sc.scalar(1e-6, unit='Å')
+
+ res = beamline.scattering_angles_with_gravity(
+ incident_beam=incident_beam,
+ scattered_beam=scattered_beam,
+ wavelength=wavelength,
+ gravity=gravity,
+ )
+
+ sc.testing.assert_allclose(
+ res['two_theta'],
+ sc.scalar(polar, unit='rad'),
+ atol=sc.scalar(1e-15, unit='rad'),
+ )
+ # We can't easily test phi here because it depends on the rotation of gravity.
+
+
+@pytest.mark.parametrize('polar', [np.pi / 3, np.pi / 2, 2 * np.pi / 3, np.pi])
+@pytest.mark.parametrize('azimuthal', [0.0, np.pi / 2, np.pi])
+def test_scattering_angles_with_gravity_reproduces_angles_orthogonal_coords(
+ polar: float,
+ azimuthal: float,
+):
+ # This case is unphysical but tests that the function reproduces
+ # the expected angles using a rotated vector.
+ # Here, incident_beam and gravity are orthogonal to use the optimized code.
+
gravity = sc.vector([0.0, -1e-11, 0.0], unit='cm/s^2')
incident_beam = sc.vector([0.0, 0.0, 968.0], unit='cm')
@@ -244,6 +324,42 @@ def test_scattering_angles_with_gravity_reproduces_angles_azimuth_greater_pi(
# This case is unphysical but tests that the function reproduces
# the expected angles using a rotated vector.
+ gravity = sc.vector([-2e-12, -1e-11, 1e-13], unit='cm/s^2')
+ incident_beam = sc.vector([0.0, 0.0, 968.0], unit='cm')
+
+ # With this definition, the x-axis has azimuthal=0.
+ rot1 = sc.spatial.rotations_from_rotvecs(sc.vector([-polar, 0, 0], unit='rad'))
+ rot2 = sc.spatial.rotations_from_rotvecs(
+ sc.vector([0, 0, azimuthal - np.pi / 2], unit='rad')
+ )
+ scattered_beam = rot2 * (rot1 * incident_beam)
+
+ wavelength = sc.scalar(1e-6, unit='Å')
+
+ res = beamline.scattering_angles_with_gravity(
+ incident_beam=incident_beam,
+ scattered_beam=scattered_beam,
+ wavelength=wavelength,
+ gravity=gravity,
+ )
+
+ sc.testing.assert_allclose(
+ res['two_theta'],
+ sc.scalar(polar, unit='rad'),
+ atol=sc.scalar(1e-15, unit='rad'),
+ )
+ # We can't easily test phi here because it depends on the rotation of gravity.
+
+
+@pytest.mark.parametrize('polar', [np.pi / 3, np.pi / 2, 2 * np.pi / 3, np.pi])
+@pytest.mark.parametrize('azimuthal', [3 * np.pi / 2, 8 * np.pi / 5, 2 * np.pi])
+def test_scattering_angles_with_gravity_reproduces_angles_azimuth_greater_pi_orthogonal_coords( # noqa: E501
+ polar: float, azimuthal: float
+):
+ # This case is unphysical but tests that the function reproduces
+ # the expected angles using a rotated vector.
+ # Here, incident_beam and gravity are orthogonal to use the optimized code.
+
gravity = sc.vector([0.0, -1e-11, 0.0], unit='cm/s^2')
incident_beam = sc.vector([0.0, 0.0, 968.0], unit='cm')
@@ -502,11 +618,12 @@ def test_scattering_angles_with_gravity_binned_data():
sc.testing.assert_identical(scattered_beam, original_scattered_beam)
-def test_scattering_angles_with_gravity_uses_wavelength_dtype():
+@pytest.mark.parametrize('g', [[0.0, -9.81, 0.0], [2.3, 0.4, -7.5]])
+def test_scattering_angles_with_gravity_uses_wavelength_dtype(g: list[float]):
wavelength = sc.array(
dims=['wavelength'], values=[1.6, 0.7], unit='Å', dtype='float32'
)
- gravity = sc.vector([0.0, -9.81, 0.0], unit='m/s^2')
+ gravity = sc.vector(g, unit='m/s^2')
incident_beam = sc.vector([0.0, 0.0, 41.1], unit='m')
scattered_beam = sc.vector([1.8, 2.5, 3.6], unit='m')
@@ -520,9 +637,10 @@ def test_scattering_angles_with_gravity_uses_wavelength_dtype():
assert res['phi'].dtype == 'float32'
-def test_scattering_angles_with_gravity_supports_mismatching_units():
+@pytest.mark.parametrize('g', [[0.0, -9.81, 0.0], [2.3, 0.4, -7.5]])
+def test_scattering_angles_with_gravity_supports_mismatching_units(g: list[float]):
wavelength = sc.array(dims=['wavelength'], values=[1.6, 0.7], unit='Å')
- gravity = sc.vector([0.0, -9810, 0.0], unit='m/ms^2')
+ gravity = sc.vector(g, unit='m/ms^2')
incident_beam = sc.vector([0.0, 0.0, 410], unit='cm')
scattered_beam = sc.vector([180, 1800, 2400], unit='mm')
@@ -533,17 +651,50 @@ def test_scattering_angles_with_gravity_supports_mismatching_units():
gravity=gravity,
)
- expected = _reference_scattering_angles_with_gravity(
+ expected = beamline.scattering_angles_with_gravity(
incident_beam=incident_beam.to(unit='m'),
scattered_beam=scattered_beam.to(unit='m'),
wavelength=wavelength,
gravity=gravity.to(unit='m/s^2'),
)
- sc.testing.assert_allclose(res['two_theta'], expected['two_theta'])
+ sc.testing.assert_allclose(
+ res['two_theta'], expected['two_theta'], rtol=sc.scalar(1e-6)
+ )
sc.testing.assert_allclose(res['phi'], expected['phi'])
+@given(rotation=rotations())
+@settings(max_examples=10)
+def test_scattering_angles_with_gravity_invariant_under_rotation(
+ rotation: sc.Variable,
+) -> None:
+ wavelength = sc.array(dims=['wavelength'], values=[1.6, 0.9, 0.7], unit='Å')
+ # Gravity and incident_beam are not aligned with the coordinate system
+ # but orthogonal to each other.
+ gravity = sc.vector([-0.3, -9.81, 0.01167883211678832], unit='m/s^2')
+ incident_beam = sc.vector([1.6, 0.0, 41.1], unit='m')
+ scattered_beam = sc.vectors(
+ dims=['det'], values=[[1.8, 2.5, 3.6], [-0.4, -1.7, 2.9]], unit='m'
+ )
+ original = beamline.scattering_angles_with_gravity(
+ incident_beam=incident_beam,
+ scattered_beam=scattered_beam,
+ wavelength=wavelength,
+ gravity=gravity,
+ )
+
+ rotated = beamline.scattering_angles_with_gravity(
+ incident_beam=rotation * incident_beam,
+ scattered_beam=rotation * scattered_beam,
+ wavelength=wavelength,
+ gravity=rotation * gravity,
+ )
+
+ sc.testing.assert_allclose(original['two_theta'], rotated['two_theta'])
+ sc.testing.assert_allclose(original['phi'], rotated['phi'])
+
+
def test_scattering_angle_in_yz_plane_requires_gravity_orthogonal_to_incident_beam():
incident_beam = sc.vector([0.564, 1.2, -10.4], unit='m')
scattered_beam = sc.vectors(
@@ -890,13 +1041,3 @@ def test_beam_aligned_unit_vectors_complicated_inputs():
sc.testing.assert_allclose(sc.dot(ex, ey), sc.scalar(0.0), atol=sc.scalar(1e-16))
sc.testing.assert_allclose(sc.dot(ey, ez), sc.scalar(0.0), atol=sc.scalar(1e-16))
sc.testing.assert_allclose(sc.dot(ez, ex), sc.scalar(0.0), atol=sc.scalar(1e-16))
-
-
-def test_beam_aligned_unit_vectors_requires_orthogonal_inputs():
- with pytest.raises(
- ValueError, match='`gravity` and `incident_beam` must be orthogonal'
- ):
- beamline.beam_aligned_unit_vectors(
- incident_beam=sc.vector([0.0, 0.0, 3.1], unit='mm'),
- gravity=sc.vector([0.0, -4.6, 1.0], unit='m/s/s'),
- )