From a18279f58547b7e9dfd21b07c340f31273ca2fac Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Fri, 9 Aug 2024 09:36:20 +0200 Subject: [PATCH 01/12] Compute unit vectors for non-orthogonal inputs --- src/scippneutron/conversion/beamline.py | 20 ++++----- tests/conversion/beamline_conversions_test.py | 42 +++++++++++++++---- 2 files changed, 41 insertions(+), 21 deletions(-) diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index 65019e3e0..b1ab2b486 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -71,7 +71,7 @@ .. math:: - \hat{e}_z &= b_1 / |b_1| \\ + \hat{e}_z &= b_1 - (b_1 \cdot \hat{e}_y) \hat{e}_y \\ \hat{e}_y &= -g / |g| \\ \hat{e}_x &= \hat{e}_y \times \hat{e}_z @@ -349,7 +349,7 @@ def beam_aligned_unit_vectors( .. math:: - \hat{e}_z &= b_1 / |b_1| \\ + \hat{e}_z &= b_1 - (b_1 \cdot \hat{e}_y) \hat{e}_y \\ \hat{e}_y &= -g / |g| \\ \hat{e}_x &= \hat{e}_y \times \hat{e}_z @@ -372,18 +372,12 @@ def beam_aligned_unit_vectors( 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) - ): - 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.' - ) - - ez = incident_beam / sc.norm(incident_beam) ey = -gravity / sc.norm(gravity) + + # project incoming_beam onto the plane perpendicular to gravity + z = incident_beam - sc.dot(incident_beam, ey) * ey + ez = z / sc.norm(z) + ex = sc.cross(ey, ez) return { 'beam_aligned_unit_x': ex, diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py index 0039974ae..8369af9a0 100644 --- a/tests/conversion/beamline_conversions_test.py +++ b/tests/conversion/beamline_conversions_test.py @@ -892,11 +892,37 @@ def test_beam_aligned_unit_vectors_complicated_inputs(): 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'), - ) +def test_beam_aligned_unit_vectors_simple_non_orthogonal_inputs(): + incident_beam = sc.vector([0.0, -1.1, 14.7], unit='m') + gravity = sc.vector([0.0, -9.3, 0.0], unit='m/s/s') + res = beamline.beam_aligned_unit_vectors( + incident_beam=incident_beam, gravity=gravity + ) + sc.testing.assert_allclose(res['beam_aligned_unit_x'], sc.vector([1.0, 0.0, 0.0])) + sc.testing.assert_allclose(res['beam_aligned_unit_y'], sc.vector([0.0, 1.0, 0.0])) + sc.testing.assert_allclose(res['beam_aligned_unit_z'], sc.vector([0.0, 0.0, 1.0])) + + +def test_beam_aligned_unit_vectors_complex_non_orthogonal_inputs(): + incident_beam = sc.vector([2.3, -1.1, 14.7], unit='m') + gravity = sc.vector([0.0, -9.3, 0.0], unit='m/s/s') + res = beamline.beam_aligned_unit_vectors( + incident_beam=incident_beam, gravity=gravity + ) + z = sc.vector([2.3, 0.0, 14.7], unit='m') + ez = z / sc.norm(z) + sc.testing.assert_allclose(sc.dot(res['beam_aligned_unit_x'], ez), sc.scalar(0.0)) + sc.testing.assert_allclose(sc.dot(res['beam_aligned_unit_y'], ez), sc.scalar(0.0)) + sc.testing.assert_allclose(res['beam_aligned_unit_y'], sc.vector([0.0, 1.0, 0.0])) + sc.testing.assert_allclose(res['beam_aligned_unit_z'], ez) + + +def test_beam_aligned_unit_vectors_non_orthogonal_inputs_gravity_along_z(): + incident_beam = sc.vector([0.0, 6.7, 0.3], unit='m') + gravity = sc.vector([0.0, 0.0, 3.4], unit='m/s/s') + res = beamline.beam_aligned_unit_vectors( + incident_beam=incident_beam, gravity=gravity + ) + sc.testing.assert_allclose(res['beam_aligned_unit_x'], sc.vector([1.0, 0.0, 0.0])) + sc.testing.assert_allclose(res['beam_aligned_unit_y'], sc.vector([0.0, 0.0, -1.0])) + sc.testing.assert_allclose(res['beam_aligned_unit_z'], sc.vector([0.0, 1.0, 0.0])) From 4982020e9460fabf5d0c4319b73c02b443ecaab5 Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Fri, 9 Aug 2024 10:38:53 +0200 Subject: [PATCH 02/12] Test two_theta invariance under rotation --- tests/conversion/beamline_conversions_test.py | 35 +++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py index 8369af9a0..0ded8911b 100644 --- a/tests/conversion/beamline_conversions_test.py +++ b/tests/conversion/beamline_conversions_test.py @@ -3,10 +3,14 @@ # @author Jan-Lukas Wynen import numpy as np +import numpy.typing as npt import pytest 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 @@ -157,6 +161,37 @@ def test_two_theta_invariant_under_reflection_about_incident_beam(): sc.testing.assert_allclose(two_theta[0], two_theta[1]) +def quaternions() -> st.SearchStrategy[npt.NDArray[float]]: + 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: q.all()) + .map(lambda q: q / np.linalg.norm(q)) + ) + + +@given(q=quaternions()) +@settings(max_examples=10) +def test_two_theta_invariant_under_rotation(q: npt.NDArray[float]) -> 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 + ) + + rotation = sc.spatial.rotation(value=q) + rotated = beamline.two_theta( + incident_beam=rotation * incident_beam, scattered_beam=rotation * scattered_beam + ) + + sc.testing.assert_allclose(original, rotated) + + +# TODO def test_scattering_angles_requires_gravity_orthogonal_to_incident_beam(): incident_beam = sc.vector([0.564, 1.2, -10.4], unit='m') scattered_beam = sc.vectors( From 39b00b32d7aee0066280d6301a048b5de10c69e7 Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Fri, 9 Aug 2024 10:47:06 +0200 Subject: [PATCH 03/12] Test angles with gravity invariance under rotation --- tests/conversion/beamline_conversions_test.py | 60 ++++++++++++++----- 1 file changed, 45 insertions(+), 15 deletions(-) diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py index 0ded8911b..ce248dc88 100644 --- a/tests/conversion/beamline_conversions_test.py +++ b/tests/conversion/beamline_conversions_test.py @@ -3,7 +3,6 @@ # @author Jan-Lukas Wynen import numpy as np -import numpy.typing as npt import pytest import scipp as sc import scipp.constants @@ -15,6 +14,17 @@ 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' @@ -161,20 +171,9 @@ def test_two_theta_invariant_under_reflection_about_incident_beam(): sc.testing.assert_allclose(two_theta[0], two_theta[1]) -def quaternions() -> st.SearchStrategy[npt.NDArray[float]]: - 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: q.all()) - .map(lambda q: q / np.linalg.norm(q)) - ) - - -@given(q=quaternions()) +@given(rotation=rotations()) @settings(max_examples=10) -def test_two_theta_invariant_under_rotation(q: npt.NDArray[float]) -> None: +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' @@ -183,7 +182,6 @@ def test_two_theta_invariant_under_rotation(q: npt.NDArray[float]) -> None: incident_beam=incident_beam, scattered_beam=scattered_beam ) - rotation = sc.spatial.rotation(value=q) rotated = beamline.two_theta( incident_beam=rotation * incident_beam, scattered_beam=rotation * scattered_beam ) @@ -579,6 +577,38 @@ def test_scattering_angles_with_gravity_supports_mismatching_units(): 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']) + + +# TODO 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( From 3f4a6baabb2c39885956a729599364109b9bd101 Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Fri, 23 Aug 2024 14:48:15 +0200 Subject: [PATCH 04/12] Use fallback impl instead of trying to fix old one --- src/scippneutron/conversion/beamline.py | 79 +++++++- tests/conversion/beamline_conversions_test.py | 168 ++++++++++++------ 2 files changed, 189 insertions(+), 58 deletions(-) diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index b1ab2b486..ddc1e54b9 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -44,6 +44,9 @@ However, when we need ``phi``, or need to correct ``two_theta`` for gravity, we need the actual coordinate system. +TODO change def of coord system +TODO check whether angle computation is correct: + should be between rays, not plane ScippNeutron uses a coordinate system aligned with the incident beam and gravity. ScippNeutron's coordinate system corresponds to that of @@ -349,7 +352,7 @@ def beam_aligned_unit_vectors( .. math:: - \hat{e}_z &= b_1 - (b_1 \cdot \hat{e}_y) \hat{e}_y \\ + \hat{e}_z &= b_1 / |b_1| \\ \hat{e}_y &= -g / |g| \\ \hat{e}_x &= \hat{e}_y \times \hat{e}_z @@ -372,12 +375,18 @@ def beam_aligned_unit_vectors( 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) + ): + 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.' + ) + + ez = incident_beam / sc.norm(incident_beam) ey = -gravity / sc.norm(gravity) - - # project incoming_beam onto the plane perpendicular to gravity - z = incident_beam - sc.dot(incident_beam, ey) * ey - ez = z / sc.norm(z) - ex = sc.cross(ey, ez) return { 'beam_aligned_unit_x': ex, @@ -433,6 +442,7 @@ class SphericalCoordinates(TypedDict): """ +# TODO docs: show simple formula and explain complicated one for optimisation def scattering_angles_with_gravity( incident_beam: sc.Variable, scattered_beam: sc.Variable, @@ -511,9 +521,62 @@ def scattering_angles_with_gravity( Ignores the ``x`` component when computing ``theta``. This is used in reflectometry. """ - unit_vectors = beam_aligned_unit_vectors( - incident_beam=incident_beam, gravity=gravity + try: + unit_vectors = beam_aligned_unit_vectors( + incident_beam=incident_beam, gravity=gravity + ) + except ValueError as err: + if "orthogonal" not in err.args[0]: + raise + 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( + unit_vectors=unit_vectors, + 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: + # TODO check works for any units + # TODO check uses dtype of wavelength + drop_distance = _drop_due_to_gravity( + distance=sc.norm(scattered_beam), wavelength=wavelength, gravity=gravity + ) + drop = drop_distance * (gravity / sc.norm(gravity)) + drop += scattered_beam + return { + 'two_theta': two_theta(incident_beam=incident_beam, scattered_beam=drop), + 'phi': None, + # TODO: how to define phi? + # it is define wrt. a coord system, so we need one + # we could use the idea of projecting `incident_beam` onto a plane + # that is perpendicular to `gravity`. This defines a system that should + # be close enough to the lab system of the instrument: + # # project incoming_beam onto the plane perpendicular to gravity + # z = incident_beam - sc.dot(incident_beam, ey) * ey + # ez = z / sc.norm(z) + } + + +# 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( + unit_vectors: BeamAlignedUnitVectors, + scattered_beam: sc.Variable, + wavelength: sc.Variable, + gravity: sc.Variable, +) -> SphericalCoordinates: ex = unit_vectors['beam_aligned_unit_x'] ey = unit_vectors['beam_aligned_unit_y'] ez = unit_vectors['beam_aligned_unit_z'] diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py index ce248dc88..77c8428da 100644 --- a/tests/conversion/beamline_conversions_test.py +++ b/tests/conversion/beamline_conversions_test.py @@ -189,28 +189,38 @@ def test_two_theta_invariant_under_rotation(rotation: sc.Variable) -> None: sc.testing.assert_allclose(original, rotated) -# TODO -def test_scattering_angles_requires_gravity_orthogonal_to_incident_beam(): - incident_beam = sc.vector([0.564, 1.2, -10.4], unit='m') +@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]], unit='m' + 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'], @@ -234,12 +244,53 @@ def test_scattering_angles_with_gravity_small_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]) +@given(gravity_rotation=rotations()) +@settings(max_examples=20) def test_scattering_angles_with_gravity_reproduces_angles( - polar: float, azimuthal: float + polar: float, azimuthal: float, gravity_rotation: sc.Variable ): # This case is unphysical but tests that the function reproduces # the expected angles using a rotated vector. + gravity = gravity_rotation * sc.vector([0.0, -1e-11, 0.0], 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'), + ) + sc.testing.assert_allclose( + res['phi'], sc.scalar(azimuthal, unit='rad'), atol=sc.scalar(1e-15, unit='rad') + ) + + +@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') @@ -271,11 +322,54 @@ def test_scattering_angles_with_gravity_reproduces_angles( @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]) +@given(gravity_rotation=rotations()) +@settings(max_examples=20) def test_scattering_angles_with_gravity_reproduces_angles_azimuth_greater_pi( + polar: float, azimuthal: float, gravity_rotation: sc.Variable +): + # This case is unphysical but tests that the function reproduces + # the expected angles using a rotated vector. + + gravity = gravity_rotation * sc.vector([0.0, -1e-11, 0.0], 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'), + ) + # phi has range (-pi, pi], so for azimuthal > pi, we get a negative result + sc.testing.assert_allclose( + res['phi'], + sc.scalar(azimuthal - 2 * np.pi, unit='rad'), + atol=sc.scalar(1e-15, unit='rad'), + ) + + +@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') @@ -957,37 +1051,11 @@ def test_beam_aligned_unit_vectors_complicated_inputs(): sc.testing.assert_allclose(sc.dot(ez, ex), sc.scalar(0.0), atol=sc.scalar(1e-16)) -def test_beam_aligned_unit_vectors_simple_non_orthogonal_inputs(): - incident_beam = sc.vector([0.0, -1.1, 14.7], unit='m') - gravity = sc.vector([0.0, -9.3, 0.0], unit='m/s/s') - res = beamline.beam_aligned_unit_vectors( - incident_beam=incident_beam, gravity=gravity - ) - sc.testing.assert_allclose(res['beam_aligned_unit_x'], sc.vector([1.0, 0.0, 0.0])) - sc.testing.assert_allclose(res['beam_aligned_unit_y'], sc.vector([0.0, 1.0, 0.0])) - sc.testing.assert_allclose(res['beam_aligned_unit_z'], sc.vector([0.0, 0.0, 1.0])) - - -def test_beam_aligned_unit_vectors_complex_non_orthogonal_inputs(): - incident_beam = sc.vector([2.3, -1.1, 14.7], unit='m') - gravity = sc.vector([0.0, -9.3, 0.0], unit='m/s/s') - res = beamline.beam_aligned_unit_vectors( - incident_beam=incident_beam, gravity=gravity - ) - z = sc.vector([2.3, 0.0, 14.7], unit='m') - ez = z / sc.norm(z) - sc.testing.assert_allclose(sc.dot(res['beam_aligned_unit_x'], ez), sc.scalar(0.0)) - sc.testing.assert_allclose(sc.dot(res['beam_aligned_unit_y'], ez), sc.scalar(0.0)) - sc.testing.assert_allclose(res['beam_aligned_unit_y'], sc.vector([0.0, 1.0, 0.0])) - sc.testing.assert_allclose(res['beam_aligned_unit_z'], ez) - - -def test_beam_aligned_unit_vectors_non_orthogonal_inputs_gravity_along_z(): - incident_beam = sc.vector([0.0, 6.7, 0.3], unit='m') - gravity = sc.vector([0.0, 0.0, 3.4], unit='m/s/s') - res = beamline.beam_aligned_unit_vectors( - incident_beam=incident_beam, gravity=gravity - ) - sc.testing.assert_allclose(res['beam_aligned_unit_x'], sc.vector([1.0, 0.0, 0.0])) - sc.testing.assert_allclose(res['beam_aligned_unit_y'], sc.vector([0.0, 0.0, -1.0])) - sc.testing.assert_allclose(res['beam_aligned_unit_z'], sc.vector([0.0, 1.0, 0.0])) +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'), + ) From 2e2f7e4eecb59513c475cc422dc647e572d1edbe Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Thu, 29 Aug 2024 14:08:27 +0200 Subject: [PATCH 05/12] Implement phi for non-orthogonal beams --- src/scippneutron/conversion/beamline.py | 57 ++++++++++++++----- tests/conversion/beamline_conversions_test.py | 43 +++++++++----- 2 files changed, 72 insertions(+), 28 deletions(-) diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index ddc1e54b9..a11b67fc9 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -329,7 +329,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): @@ -548,24 +554,49 @@ def _scattering_angles_with_gravity_generic( wavelength: sc.Variable, gravity: sc.Variable, ) -> SphericalCoordinates: - # TODO check works for any units - # TODO check uses dtype of wavelength + unit_vectors = _projected_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), - 'phi': None, - # TODO: how to define phi? - # it is define wrt. a coord system, so we need one - # we could use the idea of projecting `incident_beam` onto a plane - # that is perpendicular to `gravity`. This defines a system that should - # be close enough to the lab system of the instrument: - # # project incoming_beam onto the plane perpendicular to gravity - # z = incident_beam - sc.dot(incident_beam, ey) * ey - # ez = z / sc.norm(z) + 'two_theta': two_theta(incident_beam=incident_beam, scattered_beam=drop).to( + dtype=elem_dtype(wavelength), copy=False + ), + 'phi': phi, + } + + +def _projected_beam_aligned_unit_vectors( + incident_beam: sc.Variable, gravity: sc.Variable +) -> BeamAlignedUnitVectors: + 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( + "Cannot construct a coordinate system. The incident beam and " + "gravity are parallel to each other." + ) + ez = z / sc.norm(z) + ex = sc.cross(ey, ez) + return { + 'beam_aligned_unit_x': ex, + 'beam_aligned_unit_y': ey, + 'beam_aligned_unit_z': ez, } diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py index 77c8428da..4d027313e 100644 --- a/tests/conversion/beamline_conversions_test.py +++ b/tests/conversion/beamline_conversions_test.py @@ -7,7 +7,7 @@ import scipp as sc import scipp.constants import scipp.testing -from hypothesis import given, settings +from hypothesis import assume, given, settings from hypothesis import strategies as st from hypothesis.extra import numpy as npst @@ -255,6 +255,9 @@ def test_scattering_angles_with_gravity_reproduces_angles( gravity = gravity_rotation * sc.vector([0.0, -1e-11, 0.0], unit='cm/s^2') incident_beam = sc.vector([0.0, 0.0, 968.0], unit='cm') + # We cannot compute phi when gravity and incident_beam are parallel. + assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-10)) + # 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( @@ -276,9 +279,7 @@ def test_scattering_angles_with_gravity_reproduces_angles( sc.scalar(polar, unit='rad'), atol=sc.scalar(1e-15, unit='rad'), ) - sc.testing.assert_allclose( - res['phi'], sc.scalar(azimuthal, 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]) @@ -333,6 +334,9 @@ def test_scattering_angles_with_gravity_reproduces_angles_azimuth_greater_pi( gravity = gravity_rotation * sc.vector([0.0, -1e-11, 0.0], unit='cm/s^2') incident_beam = sc.vector([0.0, 0.0, 968.0], unit='cm') + # We cannot compute phi when gravity and incident_beam are parallel. + assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-10)) + # 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( @@ -354,12 +358,7 @@ def test_scattering_angles_with_gravity_reproduces_angles_azimuth_greater_pi( sc.scalar(polar, unit='rad'), atol=sc.scalar(1e-15, unit='rad'), ) - # phi has range (-pi, pi], so for azimuthal > pi, we get a negative result - sc.testing.assert_allclose( - res['phi'], - sc.scalar(azimuthal - 2 * np.pi, 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]) @@ -629,14 +628,21 @@ 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(): +@given(gravity_rotation=rotations()) +@settings(max_examples=20) +def test_scattering_angles_with_gravity_uses_wavelength_dtype( + gravity_rotation: sc.Variable, +): 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 = gravity_rotation * sc.vector([0.0, -9.81, 0.0], 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') + # We cannot compute phi when gravity and incident_beam are parallel. + assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-10)) + res = beamline.scattering_angles_with_gravity( incident_beam=incident_beam, scattered_beam=scattered_beam, @@ -647,12 +653,19 @@ def test_scattering_angles_with_gravity_uses_wavelength_dtype(): assert res['phi'].dtype == 'float32' -def test_scattering_angles_with_gravity_supports_mismatching_units(): +@given(gravity_rotation=rotations()) +@settings(max_examples=20) +def test_scattering_angles_with_gravity_supports_mismatching_units( + gravity_rotation: sc.Variable, +): wavelength = sc.array(dims=['wavelength'], values=[1.6, 0.7], unit='Å') - gravity = sc.vector([0.0, -9810, 0.0], unit='m/ms^2') + gravity = gravity_rotation * sc.vector([0.0, -9810, 0.0], unit='m/ms^2') incident_beam = sc.vector([0.0, 0.0, 410], unit='cm') scattered_beam = sc.vector([180, 1800, 2400], unit='mm') + # We cannot compute phi when gravity and incident_beam are parallel. + assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-10)) + res = beamline.scattering_angles_with_gravity( incident_beam=incident_beam, scattered_beam=scattered_beam, @@ -660,7 +673,7 @@ 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, From 37cd82387fbe47ee5e993572dd8ca73e510583d7 Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Thu, 29 Aug 2024 14:55:43 +0200 Subject: [PATCH 06/12] Update angle docs --- .../beamline/beamline_coordinates_dark.svg | 378 +++++++------- .../beamline/beamline_coordinates_light.svg | 486 +++++++++--------- src/scippneutron/conversion/beamline.py | 108 ++-- tests/conversion/beamline_conversions_test.py | 14 +- 4 files changed, 508 insertions(+), 478 deletions(-) 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 @@ xxyyppposition(detector)sample_positiongravitygravityposition(detector)sample_positionsource_positionzzxxyygravitygravityincident_beamincident_beamscattered_beamscattered_beampp + 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 a11b67fc9..b1b7994cf 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -41,19 +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``. -TODO change def of coord system -TODO check whether angle computation is correct: - should be between rays, not plane 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.) @@ -74,16 +71,15 @@ .. math:: - \hat{e}_z &= b_1 - (b_1 \cdot \hat{e}_y) \hat{e}_y \\ \hat{e}_y &= -g / |g| \\ + \hat{e}_z &= b_1 - (b_1 \cdot \hat{e}_y) \hat{e}_y \\ \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 @@ -358,8 +354,8 @@ def beam_aligned_unit_vectors( .. math:: - \hat{e}_z &= b_1 / |b_1| \\ \hat{e}_y &= -g / |g| \\ + \hat{e}_z &= b_1 - (b_1 \cdot \hat{e}_y) \hat{e}_y \\ \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. @@ -376,23 +372,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, @@ -478,18 +478,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} \\ + \mathsf{tan}(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|` @@ -527,13 +542,10 @@ def scattering_angles_with_gravity( Ignores the ``x`` component when computing ``theta``. This is used in reflectometry. """ - try: - unit_vectors = beam_aligned_unit_vectors( - incident_beam=incident_beam, gravity=gravity - ) - except ValueError as err: - if "orthogonal" not in err.args[0]: - raise + 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, @@ -541,7 +553,7 @@ def scattering_angles_with_gravity( gravity=gravity, ) return _scattering_angles_with_gravity_orthogonal_coords( - unit_vectors=unit_vectors, + incident_beam=incident_beam, scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, @@ -554,7 +566,7 @@ def _scattering_angles_with_gravity_generic( wavelength: sc.Variable, gravity: sc.Variable, ) -> SphericalCoordinates: - unit_vectors = _projected_beam_aligned_unit_vectors( + unit_vectors = beam_aligned_unit_vectors( incident_beam=incident_beam, gravity=gravity ) ex = unit_vectors['beam_aligned_unit_x'] @@ -579,35 +591,17 @@ def _scattering_angles_with_gravity_generic( } -def _projected_beam_aligned_unit_vectors( - incident_beam: sc.Variable, gravity: sc.Variable -) -> BeamAlignedUnitVectors: - 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( - "Cannot construct a coordinate system. The incident beam and " - "gravity are parallel to each other." - ) - ez = z / sc.norm(z) - ex = sc.cross(ey, ez) - return { - 'beam_aligned_unit_x': ex, - 'beam_aligned_unit_y': ey, - 'beam_aligned_unit_z': ez, - } - - # 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( - unit_vectors: BeamAlignedUnitVectors, + 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'] ez = unit_vectors['beam_aligned_unit_z'] diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py index 4d027313e..2313ba83c 100644 --- a/tests/conversion/beamline_conversions_test.py +++ b/tests/conversion/beamline_conversions_test.py @@ -680,7 +680,9 @@ def test_scattering_angles_with_gravity_supports_mismatching_units( 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']) @@ -1062,13 +1064,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'), - ) From dcc2b83018f33d2abd2c2502d444ff6c7c1a4322 Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Thu, 29 Aug 2024 15:01:57 +0200 Subject: [PATCH 07/12] Keep scattering_angle_in_yz_plane disabled for non-orthogonal inputs --- src/scippneutron/conversion/beamline.py | 10 ++++++++++ tests/conversion/beamline_conversions_test.py | 1 - 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index b1b7994cf..be31f5214 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -691,6 +691,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 2313ba83c..5fa1c46a3 100644 --- a/tests/conversion/beamline_conversions_test.py +++ b/tests/conversion/beamline_conversions_test.py @@ -717,7 +717,6 @@ def test_scattering_angles_with_gravity_invariant_under_rotation( sc.testing.assert_allclose(original['phi'], rotated['phi']) -# TODO 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( From 87e6cdb02ba6010bd16b72b52668d69f0a827f86 Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Thu, 29 Aug 2024 15:46:11 +0200 Subject: [PATCH 08/12] Use tighter constraints --- tests/conversion/beamline_conversions_test.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py index 5fa1c46a3..c194d84fc 100644 --- a/tests/conversion/beamline_conversions_test.py +++ b/tests/conversion/beamline_conversions_test.py @@ -256,7 +256,7 @@ def test_scattering_angles_with_gravity_reproduces_angles( incident_beam = sc.vector([0.0, 0.0, 968.0], unit='cm') # We cannot compute phi when gravity and incident_beam are parallel. - assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-10)) + assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-8)) # With this definition, the x-axis has azimuthal=0. rot1 = sc.spatial.rotations_from_rotvecs(sc.vector([-polar, 0, 0], unit='rad')) @@ -335,7 +335,7 @@ def test_scattering_angles_with_gravity_reproduces_angles_azimuth_greater_pi( incident_beam = sc.vector([0.0, 0.0, 968.0], unit='cm') # We cannot compute phi when gravity and incident_beam are parallel. - assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-10)) + assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-8)) # With this definition, the x-axis has azimuthal=0. rot1 = sc.spatial.rotations_from_rotvecs(sc.vector([-polar, 0, 0], unit='rad')) @@ -641,7 +641,7 @@ def test_scattering_angles_with_gravity_uses_wavelength_dtype( scattered_beam = sc.vector([1.8, 2.5, 3.6], unit='m') # We cannot compute phi when gravity and incident_beam are parallel. - assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-10)) + assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-8)) res = beamline.scattering_angles_with_gravity( incident_beam=incident_beam, @@ -664,7 +664,7 @@ def test_scattering_angles_with_gravity_supports_mismatching_units( scattered_beam = sc.vector([180, 1800, 2400], unit='mm') # We cannot compute phi when gravity and incident_beam are parallel. - assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-10)) + assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-8)) res = beamline.scattering_angles_with_gravity( incident_beam=incident_beam, From 8e8ef339ef8810107135a6da4ceefe59c0a78dd0 Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Thu, 29 Aug 2024 16:09:31 +0200 Subject: [PATCH 09/12] Avoid hypothesis in difficult-to-generate cases The `assume` conditions filtered out too many cases such that hypothesis complained. --- tests/conversion/beamline_conversions_test.py | 44 +++++-------------- 1 file changed, 11 insertions(+), 33 deletions(-) diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py index c194d84fc..f3d6ac087 100644 --- a/tests/conversion/beamline_conversions_test.py +++ b/tests/conversion/beamline_conversions_test.py @@ -7,7 +7,7 @@ import scipp as sc import scipp.constants import scipp.testing -from hypothesis import assume, given, settings +from hypothesis import given, settings from hypothesis import strategies as st from hypothesis.extra import numpy as npst @@ -244,20 +244,15 @@ def test_scattering_angles_with_gravity_small_gravity_orthogonal_coords() -> Non @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]) -@given(gravity_rotation=rotations()) -@settings(max_examples=20) def test_scattering_angles_with_gravity_reproduces_angles( - polar: float, azimuthal: float, gravity_rotation: sc.Variable + polar: float, azimuthal: float ): # This case is unphysical but tests that the function reproduces # the expected angles using a rotated vector. - gravity = gravity_rotation * sc.vector([0.0, -1e-11, 0.0], unit='cm/s^2') + 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') - # We cannot compute phi when gravity and incident_beam are parallel. - assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-8)) - # 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( @@ -323,20 +318,15 @@ def test_scattering_angles_with_gravity_reproduces_angles_orthogonal_coords( @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]) -@given(gravity_rotation=rotations()) -@settings(max_examples=20) def test_scattering_angles_with_gravity_reproduces_angles_azimuth_greater_pi( - polar: float, azimuthal: float, gravity_rotation: sc.Variable + polar: float, azimuthal: float ): # This case is unphysical but tests that the function reproduces # the expected angles using a rotated vector. - gravity = gravity_rotation * sc.vector([0.0, -1e-11, 0.0], unit='cm/s^2') + 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') - # We cannot compute phi when gravity and incident_beam are parallel. - assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-8)) - # 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( @@ -628,21 +618,15 @@ def test_scattering_angles_with_gravity_binned_data(): sc.testing.assert_identical(scattered_beam, original_scattered_beam) -@given(gravity_rotation=rotations()) -@settings(max_examples=20) -def test_scattering_angles_with_gravity_uses_wavelength_dtype( - gravity_rotation: sc.Variable, -): +@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 = gravity_rotation * 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') - # We cannot compute phi when gravity and incident_beam are parallel. - assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-8)) - res = beamline.scattering_angles_with_gravity( incident_beam=incident_beam, scattered_beam=scattered_beam, @@ -653,19 +637,13 @@ def test_scattering_angles_with_gravity_uses_wavelength_dtype( assert res['phi'].dtype == 'float32' -@given(gravity_rotation=rotations()) -@settings(max_examples=20) -def test_scattering_angles_with_gravity_supports_mismatching_units( - gravity_rotation: sc.Variable, -): +@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 = gravity_rotation * 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') - # We cannot compute phi when gravity and incident_beam are parallel. - assume(np.all(sc.norm(sc.cross(gravity, incident_beam)).values > 1e-8)) - res = beamline.scattering_angles_with_gravity( incident_beam=incident_beam, scattered_beam=scattered_beam, From 12e7407f45661c012e6744938fac322f2daae5df Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Tue, 3 Sep 2024 09:37:34 +0200 Subject: [PATCH 10/12] Fix 2theta definition --- src/scippneutron/conversion/beamline.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index be31f5214..430c9cdf6 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -492,7 +492,7 @@ def scattering_angles_with_gravity( .. math:: - \mathsf{tan}(2\theta) &= \sphericalangle(b_1, b_2 + \delta_y \hat{e}_y) \\ + 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 From 16e14f2e4b15ea01f5542ddff19a30062ac1e369 Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Tue, 3 Sep 2024 09:40:42 +0200 Subject: [PATCH 11/12] Fix: show that z unit vector is normalized --- src/scippneutron/conversion/beamline.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index 430c9cdf6..f63a77497 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -72,7 +72,8 @@ .. math:: \hat{e}_y &= -g / |g| \\ - \hat{e}_z &= b_1 - (b_1 \cdot \hat{e}_y) \hat{e}_y \\ + 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. @@ -355,7 +356,8 @@ def beam_aligned_unit_vectors( .. math:: \hat{e}_y &= -g / |g| \\ - \hat{e}_z &= b_1 - (b_1 \cdot \hat{e}_y) \hat{e}_y \\ + 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. From 1306f140633831c76c6796538e8a4416a3ca5b9b Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Tue, 3 Sep 2024 09:49:39 +0200 Subject: [PATCH 12/12] Remove done TODO --- src/scippneutron/conversion/beamline.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index f63a77497..bb624f1e0 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -450,7 +450,6 @@ class SphericalCoordinates(TypedDict): """ -# TODO docs: show simple formula and explain complicated one for optimisation def scattering_angles_with_gravity( incident_beam: sc.Variable, scattered_beam: sc.Variable,