From b8cb4c4023977547c8ceffcf0bcd9768e9cd3f0d Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Mon, 24 Jun 2024 13:54:47 +0200 Subject: [PATCH 1/8] Add missing prefix --- src/scippneutron/conversion/beamline.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index 367f48174..ebbdcfc1f 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -475,8 +475,8 @@ def scattering_angles_with_gravity( .. math:: - \mathsf{tan}(2\theta) &= \frac{\sqrt{x_d^2 + y_d^{\prime\, 2}}}{z} \\ - \mathsf{tan}(\phi) &= \frac{y'_d}{x} + \mathsf{tan}(2\theta) &= \frac{\sqrt{x_d^2 + y_d^{\prime\, 2}}}{z_d} \\ + \mathsf{tan}(\phi) &= \frac{y'_d}{x_d} Attention --------- From a8606932ea400ae5da68ef788e9243094946e04a Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Tue, 25 Jun 2024 09:22:37 +0200 Subject: [PATCH 2/8] Add scattering_angle_in_yz_plane --- src/scippneutron/conversion/beamline.py | 27 ++ tests/conversion/beamline_conversions_test.py | 388 ++++++++++++++++++ 2 files changed, 415 insertions(+) diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index ebbdcfc1f..a94a8f87b 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -536,3 +536,30 @@ def scattering_angles_with_gravity( two_theta_ = sc.atan2(y=y, x=z, out=y) return {'two_theta': two_theta_, 'phi': phi} + + +def scattering_angle_in_yz_plane( + incident_beam: sc.Variable, + scattered_beam: sc.Variable, + wavelength: sc.Variable, + gravity: sc.Variable, + sample_rotation: sc.Variable, +) -> sc.Variable: + match beam_aligned_unit_vectors(incident_beam=incident_beam, gravity=gravity): + case { + 'beam_aligned_unit_y': ey, + 'beam_aligned_unit_z': ez, + }: + pass + case _: + raise RuntimeError('Unexpected return value of beam_aligned_unit_vectors') + + y = _drop_due_to_gravity( + distance=sc.norm(scattered_beam), wavelength=wavelength, gravity=gravity + ) + y += sc.dot(scattered_beam, ey).to(dtype=elem_dtype(wavelength), copy=False) + y = sc.abs(y, out=y) + z = sc.dot(scattered_beam, ez).to(dtype=elem_dtype(y), copy=False) + full_angle = sc.atan2(y=y, x=z, out=y) + full_angle -= sample_rotation.to(unit=elem_unit(full_angle), copy=False) + return full_angle diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py index bbeabbdaa..fb8cda01c 100644 --- a/tests/conversion/beamline_conversions_test.py +++ b/tests/conversion/beamline_conversions_test.py @@ -544,6 +544,394 @@ def test_scattering_angles_with_gravity_supports_mismatching_units(): sc.testing.assert_allclose(res['phi'], expected['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( + dims=['beam'], values=[[13, 24, 35], [51, -42, 33]], unit='m' + ) + 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) + sample_rotation = sc.scalar(0.5, unit='rad') + + with pytest.raises( + ValueError, match='`gravity` and `incident_beam` must be orthogonal' + ): + beamline.scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + + +@pytest.mark.parametrize('omega', [0.0, -0.1, 0.4]) +def test_scattering_angle_in_yz_plane_small_gravity(omega: float): + # This case is unphysical but tests the consistency with `two_theta`. + # Note that the scattered beam must be in the x-z plane for `two_theta` + # and `scattering_angle_from_sample` to compute the same angle. + incident_beam = sc.vector([0.0, 0.0, 10.4], unit='m') + scattered_beam = sc.vectors( + dims=['beam'], + values=[[0, 24, 0], [0, -42, 0]], + unit='m', + ) + wavelength = sc.array(dims=['wavelength'], values=[1.2, 1.6, 1.8], unit='Å') + gravity = sc.vector([0, -1e-11, 0], unit=sc.constants.g.unit) + sample_rotation = sc.scalar(omega, unit='rad') + + res = beamline.scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + expected = beamline.two_theta( + incident_beam=incident_beam, scattered_beam=scattered_beam + ).broadcast(dims=['wavelength', 'beam'], shape=[3, 2]) + expected = expected - sample_rotation + sc.testing.assert_allclose(res, expected) + + +@pytest.mark.parametrize('polar', [np.pi / 3, np.pi / 2, 2 * np.pi / 3, np.pi]) +@pytest.mark.parametrize('omega', [0.0, -0.3, 0.1]) +def test_scattering_angle_in_yz_plane_reproduces_polar_angle( + polar: float, omega: float +): + # This case is unphysical but tests that the function reproduces + # the expected angles using a rotated vector. + + 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') + + # With this definition, the x-axis has azimuthal=0. + rot1 = sc.spatial.rotations_from_rotvecs(sc.vector([-polar, 0, 0], unit='rad')) + scattered_beam = rot1 * incident_beam + + wavelength = sc.scalar(1e-6, unit='Å') + sample_rotation = sc.scalar(omega, unit='rad') + + res = beamline.scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + + sc.testing.assert_allclose( + res, + sc.scalar(polar, unit='rad') - sample_rotation, + ) + + +@pytest.mark.parametrize('polar', [np.pi / 3, np.pi / 2, 2 * np.pi / 3, np.pi]) +@pytest.mark.parametrize('omega', [0.0, -0.3, 0.1]) +def test_scattering_angle_in_yz_plane_reproduces_angles_azimuth_greater_pi( + polar: float, omega: float +): + # This case is unphysical but tests that the function reproduces + # the expected angles using a rotated vector. + + 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') + + # With this definition, the x-axis has azimuthal=0. + rot1 = sc.spatial.rotations_from_rotvecs(sc.vector([-polar, 0, 0], unit='rad')) + scattered_beam = rot1 * incident_beam + + wavelength = sc.scalar(1e-6, unit='Å') + sample_rotation = sc.scalar(omega, unit='rad') + + res = beamline.scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + + sc.testing.assert_allclose( + res, + sc.scalar(polar, unit='rad') - sample_rotation, + ) + + +@pytest.mark.parametrize('polar', [np.pi / 3, np.pi / 2, 2 * np.pi / 3, np.pi]) +@pytest.mark.parametrize('omega', [0.0, -0.3, 0.1]) +@pytest.mark.parametrize('x', [0.5, 11.4, -9.7]) +def test_scattering_angle_in_yz_plane_does_not_depend_on_x( + polar: float, omega: float, x: float +): + # This case is unphysical but tests that the function reproduces + # the expected angles using a rotated vector. + + 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') + + # With this definition, the x-axis has azimuthal=0. + rot1 = sc.spatial.rotations_from_rotvecs(sc.vector([-polar, 0, 0], unit='rad')) + scattered_beam_ref = rot1 * incident_beam + scattered_beam_shift = scattered_beam_ref + sc.vector([x, 0.0, 0.0], unit='cm') + + wavelength = sc.scalar(1e-6, unit='Å') + sample_rotation = sc.scalar(omega, unit='rad') + + res_shift = beamline.scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam_shift, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + res_ref = beamline.scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam_ref, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + + sc.testing.assert_allclose(res_shift, res_ref) + + +def test_scattering_angle_in_yz_plane_drops_in_expected_direction(): + wavelength = sc.scalar(1.6, unit='Å') + gravity = sc.vector([0.0, -sc.constants.g.value, 0.0], unit=sc.constants.g.unit) + incident_beam = sc.vector([0.0, 0.0, 41.1], unit='m') + scattered_beam = sc.vectors( + dims=['det'], values=[[0.0, 2.5, 8.6], [0.0, -1.7, 6.9]], unit='m' + ) + sample_rotation = sc.scalar(0.0, unit='rad') + + with_gravity = beamline.scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + without_gravity = beamline.two_theta( + incident_beam=incident_beam, scattered_beam=scattered_beam + ) + + # The neutron was detected above the incident beam. + # So using straight paths, it looks like it scattered at a + # smaller angle (detected at a lower y) than in the real case with gravity. + assert sc.all(with_gravity[0] > without_gravity[0]).value + # The neutron was detected below the incident beam. + # So the opposite of the above comment applies. + assert sc.all(with_gravity[1] < without_gravity[1]).value + + +def test_scattering_angle_in_yz_plane_beams_aligned_with_lab_coords(): + wavelength = sc.array(dims=['wavelength'], values=[1.6, 0.9, 0.7], unit='Å') + gravity = sc.vector([0.0, -sc.constants.g.value, 0.0], unit=sc.constants.g.unit) + incident_beam = sc.vector([0.0, 0.0, 41.1], unit='m') + scattered_beam = sc.vectors( + dims=['det'], values=[[0.0, 2.5, 3.6], [0.0, -1.7, 2.9]], unit='m' + ) + sample_rotation = sc.scalar(0.3, unit='rad') + + original_wavelength = wavelength.copy() + original_gravity = gravity.copy() + original_incident_beam = incident_beam.copy() + original_scattered_beam = scattered_beam.copy() + original_sample_rotation = sample_rotation.copy() + + res = beamline.scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + + L2 = sc.norm(scattered_beam) + x = scattered_beam.fields.x + y = scattered_beam.fields.y + z = scattered_beam.fields.z + drop = ( + L2**2 + * sc.norm(gravity) + * wavelength**2 + * sc.constants.m_n**2 + / (2 * sc.constants.h**2) + ) + drop = drop.to(unit=y.unit) + true_y = y + drop + true_scattered_beam = sc.spatial.as_vectors(x, true_y, z) + expected_theta = ( + sc.asin(abs(true_y) / sc.norm(true_scattered_beam)) - sample_rotation + ) + expected_theta = expected_theta.transpose(res.dims) + + sc.testing.assert_allclose(res, expected_theta) + + sc.testing.assert_identical(wavelength, original_wavelength) + sc.testing.assert_identical(gravity, original_gravity) + sc.testing.assert_identical(incident_beam, original_incident_beam) + sc.testing.assert_identical(scattered_beam, original_scattered_beam) + sc.testing.assert_identical(sample_rotation, original_sample_rotation) + + +def _reference_scattering_angle_in_yz_plane( + incident_beam: sc.Variable, + scattered_beam: sc.Variable, + gravity: sc.Variable, + wavelength: sc.Variable, + sample_rotation: sc.Variable, +) -> sc.Variable: + # This is a simplified, independently checked implementation. + e_z = incident_beam / sc.norm(incident_beam) + e_y = -gravity / sc.norm(gravity) + + y = sc.dot(scattered_beam, e_y) + z = sc.dot(scattered_beam, e_z) + + L2 = sc.norm(scattered_beam) + drop = ( + L2**2 + * wavelength**2 + * sc.norm(gravity) + * (sc.constants.m_n**2 / (2 * sc.constants.h**2)) + ) + dropped_y = y + drop.to(unit=y.unit) + + return sc.atan2(y=abs(dropped_y), x=z) - sample_rotation + + +def test_scattering_angle_in_yz_plane_beams_unaligned_with_lab_coords(): + 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' + ) + sample_rotation = sc.scalar(0.4, unit='rad') + + original_wavelength = wavelength.copy() + original_gravity = gravity.copy() + original_incident_beam = incident_beam.copy() + original_scattered_beam = scattered_beam.copy() + original_sample_rotation = sample_rotation.copy() + + res = beamline.scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + expected = _reference_scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + + sc.testing.assert_allclose(res, expected.transpose(res.dims)) + + sc.testing.assert_identical(wavelength, original_wavelength) + sc.testing.assert_identical(gravity, original_gravity) + sc.testing.assert_identical(incident_beam, original_incident_beam) + sc.testing.assert_identical(scattered_beam, original_scattered_beam) + sc.testing.assert_identical(sample_rotation, original_sample_rotation) + + +def test_scattering_angle_in_yz_plane_binned_data(): + wavelength = sc.array(dims=['wavelength'], values=[1.6, 0.9, 0.7], unit='Å') + wavelength = sc.bins( + dim='wavelength', + data=wavelength, + begin=sc.array(dims=['det'], values=[0, 2], unit=None), + end=sc.array(dims=['det'], values=[2, 3], unit=None), + ) + 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' + ) + sample_rotation = sc.scalar(0.2, unit='rad') + + original_wavelength = wavelength.copy() + original_gravity = gravity.copy() + original_incident_beam = incident_beam.copy() + original_scattered_beam = scattered_beam.copy() + original_sample_rotation = sample_rotation.copy() + + res = beamline.scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + expected = _reference_scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + + sc.testing.assert_allclose(res, expected.transpose(res.dims)) + sc.testing.assert_identical(wavelength, original_wavelength) + sc.testing.assert_identical(gravity, original_gravity) + sc.testing.assert_identical(incident_beam, original_incident_beam) + sc.testing.assert_identical(scattered_beam, original_scattered_beam) + sc.testing.assert_identical(sample_rotation, original_sample_rotation) + + +def test_scattering_angle_in_yz_plane_uses_wavelength_dtype(): + 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') + incident_beam = sc.vector([0.0, 0.0, 41.1], unit='m') + scattered_beam = sc.vector([1.8, 2.5, 3.6], unit='m') + sample_rotation = sc.scalar(0.2, unit='rad') + + res = beamline.scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + assert res.dtype == 'float32' + + +def test_scattering_angle_in_yz_plane_supports_mismatching_units(): + wavelength = sc.array(dims=['wavelength'], values=[1.6, 0.7], unit='Å') + gravity = 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') + sample_rotation = sc.scalar(0.2, unit='deg') + + res = beamline.scattering_angle_in_yz_plane( + incident_beam=incident_beam, + scattered_beam=scattered_beam, + wavelength=wavelength, + gravity=gravity, + sample_rotation=sample_rotation, + ) + + expected = _reference_scattering_angle_in_yz_plane( + incident_beam=incident_beam.to(unit='m'), + scattered_beam=scattered_beam.to(unit='m'), + wavelength=wavelength, + gravity=gravity.to(unit='m/s^2'), + sample_rotation=sample_rotation.to(unit='rad'), + ) + + sc.testing.assert_allclose(res, expected) + + def test_beam_aligned_unit_vectors_axis_aligned_inputs(): res = beamline.beam_aligned_unit_vectors( incident_beam=sc.vector([0.0, 0.0, 2.1], unit='mm'), From f50d07f0605b869a5bf452a624c181f3c54f112c Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Tue, 25 Jun 2024 11:11:07 +0200 Subject: [PATCH 3/8] Remove sample_rotation from scattering_angle_in_yz_plane --- src/scippneutron/conversion/beamline.py | 5 +- tests/conversion/beamline_conversions_test.py | 64 +++---------------- 2 files changed, 9 insertions(+), 60 deletions(-) diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index a94a8f87b..f26866425 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -543,7 +543,6 @@ def scattering_angle_in_yz_plane( scattered_beam: sc.Variable, wavelength: sc.Variable, gravity: sc.Variable, - sample_rotation: sc.Variable, ) -> sc.Variable: match beam_aligned_unit_vectors(incident_beam=incident_beam, gravity=gravity): case { @@ -560,6 +559,4 @@ def scattering_angle_in_yz_plane( y += sc.dot(scattered_beam, ey).to(dtype=elem_dtype(wavelength), copy=False) y = sc.abs(y, out=y) z = sc.dot(scattered_beam, ez).to(dtype=elem_dtype(y), copy=False) - full_angle = sc.atan2(y=y, x=z, out=y) - full_angle -= sample_rotation.to(unit=elem_unit(full_angle), copy=False) - return full_angle + return sc.atan2(y=y, x=z, out=y) diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py index fb8cda01c..163549072 100644 --- a/tests/conversion/beamline_conversions_test.py +++ b/tests/conversion/beamline_conversions_test.py @@ -551,7 +551,6 @@ def test_scattering_angle_in_yz_plane_requires_gravity_orthogonal_to_incident_be ) 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) - sample_rotation = sc.scalar(0.5, unit='rad') with pytest.raises( ValueError, match='`gravity` and `incident_beam` must be orthogonal' @@ -561,12 +560,10 @@ def test_scattering_angle_in_yz_plane_requires_gravity_orthogonal_to_incident_be scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) -@pytest.mark.parametrize('omega', [0.0, -0.1, 0.4]) -def test_scattering_angle_in_yz_plane_small_gravity(omega: float): +def test_scattering_angle_in_yz_plane_small_gravity(): # This case is unphysical but tests the consistency with `two_theta`. # Note that the scattered beam must be in the x-z plane for `two_theta` # and `scattering_angle_from_sample` to compute the same angle. @@ -578,26 +575,22 @@ def test_scattering_angle_in_yz_plane_small_gravity(omega: float): ) wavelength = sc.array(dims=['wavelength'], values=[1.2, 1.6, 1.8], unit='Å') gravity = sc.vector([0, -1e-11, 0], unit=sc.constants.g.unit) - sample_rotation = sc.scalar(omega, unit='rad') res = beamline.scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) expected = beamline.two_theta( incident_beam=incident_beam, scattered_beam=scattered_beam ).broadcast(dims=['wavelength', 'beam'], shape=[3, 2]) - expected = expected - sample_rotation sc.testing.assert_allclose(res, expected) @pytest.mark.parametrize('polar', [np.pi / 3, np.pi / 2, 2 * np.pi / 3, np.pi]) -@pytest.mark.parametrize('omega', [0.0, -0.3, 0.1]) def test_scattering_angle_in_yz_plane_reproduces_polar_angle( - polar: float, omega: float + polar: float, ): # This case is unphysical but tests that the function reproduces # the expected angles using a rotated vector. @@ -610,26 +603,20 @@ def test_scattering_angle_in_yz_plane_reproduces_polar_angle( scattered_beam = rot1 * incident_beam wavelength = sc.scalar(1e-6, unit='Å') - sample_rotation = sc.scalar(omega, unit='rad') res = beamline.scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) - sc.testing.assert_allclose( - res, - sc.scalar(polar, unit='rad') - sample_rotation, - ) + sc.testing.assert_allclose(res, sc.scalar(polar, unit='rad')) @pytest.mark.parametrize('polar', [np.pi / 3, np.pi / 2, 2 * np.pi / 3, np.pi]) -@pytest.mark.parametrize('omega', [0.0, -0.3, 0.1]) def test_scattering_angle_in_yz_plane_reproduces_angles_azimuth_greater_pi( - polar: float, omega: float + polar: float, ): # This case is unphysical but tests that the function reproduces # the expected angles using a rotated vector. @@ -642,28 +629,20 @@ def test_scattering_angle_in_yz_plane_reproduces_angles_azimuth_greater_pi( scattered_beam = rot1 * incident_beam wavelength = sc.scalar(1e-6, unit='Å') - sample_rotation = sc.scalar(omega, unit='rad') res = beamline.scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) - sc.testing.assert_allclose( - res, - sc.scalar(polar, unit='rad') - sample_rotation, - ) + sc.testing.assert_allclose(res, sc.scalar(polar, unit='rad')) @pytest.mark.parametrize('polar', [np.pi / 3, np.pi / 2, 2 * np.pi / 3, np.pi]) -@pytest.mark.parametrize('omega', [0.0, -0.3, 0.1]) @pytest.mark.parametrize('x', [0.5, 11.4, -9.7]) -def test_scattering_angle_in_yz_plane_does_not_depend_on_x( - polar: float, omega: float, x: float -): +def test_scattering_angle_in_yz_plane_does_not_depend_on_x(polar: float, x: float): # This case is unphysical but tests that the function reproduces # the expected angles using a rotated vector. @@ -676,21 +655,18 @@ def test_scattering_angle_in_yz_plane_does_not_depend_on_x( scattered_beam_shift = scattered_beam_ref + sc.vector([x, 0.0, 0.0], unit='cm') wavelength = sc.scalar(1e-6, unit='Å') - sample_rotation = sc.scalar(omega, unit='rad') res_shift = beamline.scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam_shift, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) res_ref = beamline.scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam_ref, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) sc.testing.assert_allclose(res_shift, res_ref) @@ -703,14 +679,12 @@ def test_scattering_angle_in_yz_plane_drops_in_expected_direction(): scattered_beam = sc.vectors( dims=['det'], values=[[0.0, 2.5, 8.6], [0.0, -1.7, 6.9]], unit='m' ) - sample_rotation = sc.scalar(0.0, unit='rad') with_gravity = beamline.scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) without_gravity = beamline.two_theta( incident_beam=incident_beam, scattered_beam=scattered_beam @@ -732,20 +706,17 @@ def test_scattering_angle_in_yz_plane_beams_aligned_with_lab_coords(): scattered_beam = sc.vectors( dims=['det'], values=[[0.0, 2.5, 3.6], [0.0, -1.7, 2.9]], unit='m' ) - sample_rotation = sc.scalar(0.3, unit='rad') original_wavelength = wavelength.copy() original_gravity = gravity.copy() original_incident_beam = incident_beam.copy() original_scattered_beam = scattered_beam.copy() - original_sample_rotation = sample_rotation.copy() res = beamline.scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) L2 = sc.norm(scattered_beam) @@ -762,9 +733,7 @@ def test_scattering_angle_in_yz_plane_beams_aligned_with_lab_coords(): drop = drop.to(unit=y.unit) true_y = y + drop true_scattered_beam = sc.spatial.as_vectors(x, true_y, z) - expected_theta = ( - sc.asin(abs(true_y) / sc.norm(true_scattered_beam)) - sample_rotation - ) + expected_theta = sc.asin(abs(true_y) / sc.norm(true_scattered_beam)) expected_theta = expected_theta.transpose(res.dims) sc.testing.assert_allclose(res, expected_theta) @@ -773,7 +742,6 @@ def test_scattering_angle_in_yz_plane_beams_aligned_with_lab_coords(): sc.testing.assert_identical(gravity, original_gravity) sc.testing.assert_identical(incident_beam, original_incident_beam) sc.testing.assert_identical(scattered_beam, original_scattered_beam) - sc.testing.assert_identical(sample_rotation, original_sample_rotation) def _reference_scattering_angle_in_yz_plane( @@ -781,7 +749,6 @@ def _reference_scattering_angle_in_yz_plane( scattered_beam: sc.Variable, gravity: sc.Variable, wavelength: sc.Variable, - sample_rotation: sc.Variable, ) -> sc.Variable: # This is a simplified, independently checked implementation. e_z = incident_beam / sc.norm(incident_beam) @@ -799,7 +766,7 @@ def _reference_scattering_angle_in_yz_plane( ) dropped_y = y + drop.to(unit=y.unit) - return sc.atan2(y=abs(dropped_y), x=z) - sample_rotation + return sc.atan2(y=abs(dropped_y), x=z) def test_scattering_angle_in_yz_plane_beams_unaligned_with_lab_coords(): @@ -811,27 +778,23 @@ def test_scattering_angle_in_yz_plane_beams_unaligned_with_lab_coords(): scattered_beam = sc.vectors( dims=['det'], values=[[1.8, 2.5, 3.6], [-0.4, -1.7, 2.9]], unit='m' ) - sample_rotation = sc.scalar(0.4, unit='rad') original_wavelength = wavelength.copy() original_gravity = gravity.copy() original_incident_beam = incident_beam.copy() original_scattered_beam = scattered_beam.copy() - original_sample_rotation = sample_rotation.copy() res = beamline.scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) expected = _reference_scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) sc.testing.assert_allclose(res, expected.transpose(res.dims)) @@ -840,7 +803,6 @@ def test_scattering_angle_in_yz_plane_beams_unaligned_with_lab_coords(): sc.testing.assert_identical(gravity, original_gravity) sc.testing.assert_identical(incident_beam, original_incident_beam) sc.testing.assert_identical(scattered_beam, original_scattered_beam) - sc.testing.assert_identical(sample_rotation, original_sample_rotation) def test_scattering_angle_in_yz_plane_binned_data(): @@ -856,27 +818,23 @@ def test_scattering_angle_in_yz_plane_binned_data(): scattered_beam = sc.vectors( dims=['det'], values=[[1.8, 2.5, 3.6], [-0.4, -1.7, 2.9]], unit='m' ) - sample_rotation = sc.scalar(0.2, unit='rad') original_wavelength = wavelength.copy() original_gravity = gravity.copy() original_incident_beam = incident_beam.copy() original_scattered_beam = scattered_beam.copy() - original_sample_rotation = sample_rotation.copy() res = beamline.scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) expected = _reference_scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) sc.testing.assert_allclose(res, expected.transpose(res.dims)) @@ -884,7 +842,6 @@ def test_scattering_angle_in_yz_plane_binned_data(): sc.testing.assert_identical(gravity, original_gravity) sc.testing.assert_identical(incident_beam, original_incident_beam) sc.testing.assert_identical(scattered_beam, original_scattered_beam) - sc.testing.assert_identical(sample_rotation, original_sample_rotation) def test_scattering_angle_in_yz_plane_uses_wavelength_dtype(): @@ -894,14 +851,12 @@ def test_scattering_angle_in_yz_plane_uses_wavelength_dtype(): gravity = 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') - sample_rotation = sc.scalar(0.2, unit='rad') res = beamline.scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) assert res.dtype == 'float32' @@ -911,14 +866,12 @@ def test_scattering_angle_in_yz_plane_supports_mismatching_units(): gravity = 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') - sample_rotation = sc.scalar(0.2, unit='deg') res = beamline.scattering_angle_in_yz_plane( incident_beam=incident_beam, scattered_beam=scattered_beam, wavelength=wavelength, gravity=gravity, - sample_rotation=sample_rotation, ) expected = _reference_scattering_angle_in_yz_plane( @@ -926,7 +879,6 @@ def test_scattering_angle_in_yz_plane_supports_mismatching_units(): scattered_beam=scattered_beam.to(unit='m'), wavelength=wavelength, gravity=gravity.to(unit='m/s^2'), - sample_rotation=sample_rotation.to(unit='rad'), ) sc.testing.assert_allclose(res, expected) From 1692915579bd12911aae4835c33012c8293fdbbe Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Tue, 25 Jun 2024 11:13:32 +0200 Subject: [PATCH 4/8] Simplify --- src/scippneutron/conversion/beamline.py | 28 ++++++++++--------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index f26866425..b7e612e47 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -508,15 +508,12 @@ def scattering_angles_with_gravity( A dict containing the polar scattering angle ``'two_theta'`` and the azimuthal angle ``'phi'``. """ - match beam_aligned_unit_vectors(incident_beam=incident_beam, gravity=gravity): - case { - 'beam_aligned_unit_x': ex, - 'beam_aligned_unit_y': ey, - 'beam_aligned_unit_z': ez, - }: - pass - case _: - raise RuntimeError('Unexpected return value of beam_aligned_unit_vectors') + 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'] y = _drop_due_to_gravity( distance=sc.norm(scattered_beam), wavelength=wavelength, gravity=gravity @@ -544,14 +541,11 @@ def scattering_angle_in_yz_plane( wavelength: sc.Variable, gravity: sc.Variable, ) -> sc.Variable: - match beam_aligned_unit_vectors(incident_beam=incident_beam, gravity=gravity): - case { - 'beam_aligned_unit_y': ey, - 'beam_aligned_unit_z': ez, - }: - pass - case _: - raise RuntimeError('Unexpected return value of beam_aligned_unit_vectors') + unit_vectors = beam_aligned_unit_vectors( + incident_beam=incident_beam, gravity=gravity + ) + ey = unit_vectors['beam_aligned_unit_y'] + ez = unit_vectors['beam_aligned_unit_z'] y = _drop_due_to_gravity( distance=sc.norm(scattered_beam), wavelength=wavelength, gravity=gravity From 3a4ce8a3a1f814e638ffb62aad0694a61b567860 Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Tue, 25 Jun 2024 11:26:55 +0200 Subject: [PATCH 5/8] Document scattering_angle_in_yz_plane --- src/scippneutron/conversion/beamline.py | 80 ++++++++++++++++++++++--- 1 file changed, 71 insertions(+), 9 deletions(-) diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index b7e612e47..ad82fa7b8 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -480,16 +480,16 @@ def scattering_angles_with_gravity( Attention --------- - The above equation for :math:`y'_d` contains :math:`L_2^{\prime\, 2} = |b'_2|` - which in turn depends on :math:`y'_d`. - Solving this equation for :math:`y'_d` is too difficult. - Instead, we approximate :math:`L'_2 \approx L_2`. - The impact of this approximation on :math:`2\theta` is of the order of - :math:`10^{-6}` or less for beamlines at ESS. - This is within the expected statistical uncertainties and can be ignored. + The above equation for :math:`y'_d` contains :math:`L_2^{\prime\, 2} = |b'_2|` + which in turn depends on :math:`y'_d`. + Solving this equation for :math:`y'_d` is too difficult. + Instead, we approximate :math:`L'_2 \approx L_2`. + The impact of this approximation on :math:`2\theta` is of the order of + :math:`10^{-6}` or less for beamlines at ESS. + This is within the expected statistical uncertainties and can be ignored. - See the file ``docs/user-guide/auxiliary/gravity_approx/gravity_approx.typ`` - for details. + See the file ``docs/user-guide/auxiliary/gravity_approx/gravity_approx.typ`` + for details. Parameters ---------- @@ -507,6 +507,12 @@ def scattering_angles_with_gravity( : A dict containing the polar scattering angle ``'two_theta'`` and the azimuthal angle ``'phi'``. + + See also + -------- + scattering_angle_in_yz_plane: + Ignores the ``x`` component when computing ``theta``. + This is used in reflectometry. """ unit_vectors = beam_aligned_unit_vectors( incident_beam=incident_beam, gravity=gravity @@ -541,6 +547,62 @@ def scattering_angle_in_yz_plane( wavelength: sc.Variable, gravity: sc.Variable, ) -> sc.Variable: + r"""Compute polar scattering angles in the y-z plane using gravity. + + Note + ---- + This function uses the reflectometry definition of the polar scattering angle. + Other techniques define the angle w.r.t. the incident beam. + See :func:`scattering_angles_with_gravity` for those use cases. + + With the definitions given in :func:`scattering_angles_with_gravity`, + and ignoring :math:`x_d`, we get + + .. math:: + + \mathsf{tan}(\gamma) = \frac{|y_d^{\prime}|}{z_d} + + with + + .. math:: + + y'_d = y_d + \frac{|g| m_n^2}{2 h^2} L_2^{\prime\, 2} \lambda^2 + + Attention + --------- + The above equation for :math:`y'_d` contains :math:`L_2^{\prime\, 2} = |b'_2|` + which in turn depends on :math:`y'_d`. + Solving this equation for :math:`y'_d` is too difficult. + Instead, we approximate :math:`L'_2 \approx L_2`. + The impact of this approximation on :math:`\gamma` is of the order of + :math:`10^{-6}` or less for beamlines at ESS. + This is within the expected statistical uncertainties and can be ignored. + + See the file ``docs/user-guide/auxiliary/gravity_approx/gravity_approx.typ`` + for details. + + Parameters + ---------- + incident_beam: + Beam from source to sample. Expects ``dtype=vector3``. + scattered_beam: + Beam from sample to detector. Expects ``dtype=vector3``. + wavelength: + Wavelength of neutrons. + gravity: + Gravity vector. + + Returns + ------- + : + The polar scattering angle :math:`\gamma`. + + See also + -------- + scattering_angles_with_gravity: + Includes the ``x`` component when computing ``theta``. + This is used in techniques other than reflectometry. + """ unit_vectors = beam_aligned_unit_vectors( incident_beam=incident_beam, gravity=gravity ) From 3c6bc04fc647eee72b9ff2a4fb51ba0e673340f4 Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Tue, 25 Jun 2024 11:32:24 +0200 Subject: [PATCH 6/8] Link to correct notebook --- 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 ad82fa7b8..65d001e1d 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -488,7 +488,8 @@ def scattering_angles_with_gravity( :math:`10^{-6}` or less for beamlines at ESS. This is within the expected statistical uncertainties and can be ignored. - See the file ``docs/user-guide/auxiliary/gravity_approx/gravity_approx.typ`` + See `two_theta gravity correction + <../../user-guide/algorithms-background/two_theta-gravity-correction.rst>`_ for details. Parameters @@ -578,7 +579,8 @@ def scattering_angle_in_yz_plane( :math:`10^{-6}` or less for beamlines at ESS. This is within the expected statistical uncertainties and can be ignored. - See the file ``docs/user-guide/auxiliary/gravity_approx/gravity_approx.typ`` + See `two_theta gravity correction + <../../user-guide/algorithms-background/two_theta-gravity-correction.rst>`_ for details. Parameters From 976a7da1739354ace323b00ca38d3f53d982a665 Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Tue, 25 Jun 2024 11:38:08 +0200 Subject: [PATCH 7/8] Cite reflectometry paper --- docs/bibliography.bib | 13 +++++++++++++ src/scippneutron/conversion/beamline.py | 4 ++++ 2 files changed, 17 insertions(+) diff --git a/docs/bibliography.bib b/docs/bibliography.bib index bb8e39886..7211107de 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -48,3 +48,16 @@ @misc{kahan:2000:CPR note = "Lecture notes for Math 128.", URL = "http://www.cs.berkeley.edu/~wkahan/Math128/Cross.pdf", } + +@article{STAHN201644, + title = {Focusing neutron reflectometry: Implementation and experience on the TOF-reflectometer Amor}, + journal = {Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment}, + volume = {821}, + pages = {44-54}, + year = {2016}, + issn = {0168-9002}, + doi = {https://doi.org/10.1016/j.nima.2016.03.007}, + url = {https://www.sciencedirect.com/science/article/pii/S0168900216300250}, + author = {J. Stahn and A. Glavic}, + keywords = {Reflectometry, Neutron, Focusing}, +} diff --git a/src/scippneutron/conversion/beamline.py b/src/scippneutron/conversion/beamline.py index 65d001e1d..65019e3e0 100644 --- a/src/scippneutron/conversion/beamline.py +++ b/src/scippneutron/conversion/beamline.py @@ -104,6 +104,8 @@ These definitions assume that gravity can be neglected. See :func:`scattering_angles_with_gravity` for definitions that account for gravity. +And :func:`scattering_angle_in_yz_plane` for the definition used in reflectometry --- +which also includes gravity. """ from typing import TypedDict @@ -569,6 +571,8 @@ def scattering_angle_in_yz_plane( y'_d = y_d + \frac{|g| m_n^2}{2 h^2} L_2^{\prime\, 2} \lambda^2 + The angle :math:`\gamma` is defined as in Fig. 5 of :cite:`STAHN201644`. + Attention --------- The above equation for :math:`y'_d` contains :math:`L_2^{\prime\, 2} = |b'_2|` From dfcd03551f5964bd9378774830dc23aef4bcff2e Mon Sep 17 00:00:00 2001 From: Jan-Lukas Wynen Date: Wed, 26 Jun 2024 09:34:14 +0200 Subject: [PATCH 8/8] Remove redundant test --- tests/conversion/beamline_conversions_test.py | 26 ------------------- 1 file changed, 26 deletions(-) diff --git a/tests/conversion/beamline_conversions_test.py b/tests/conversion/beamline_conversions_test.py index 163549072..0039974ae 100644 --- a/tests/conversion/beamline_conversions_test.py +++ b/tests/conversion/beamline_conversions_test.py @@ -614,32 +614,6 @@ def test_scattering_angle_in_yz_plane_reproduces_polar_angle( sc.testing.assert_allclose(res, sc.scalar(polar, unit='rad')) -@pytest.mark.parametrize('polar', [np.pi / 3, np.pi / 2, 2 * np.pi / 3, np.pi]) -def test_scattering_angle_in_yz_plane_reproduces_angles_azimuth_greater_pi( - polar: float, -): - # This case is unphysical but tests that the function reproduces - # the expected angles using a rotated vector. - - 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') - - # With this definition, the x-axis has azimuthal=0. - rot1 = sc.spatial.rotations_from_rotvecs(sc.vector([-polar, 0, 0], unit='rad')) - scattered_beam = rot1 * incident_beam - - wavelength = sc.scalar(1e-6, unit='Å') - - res = beamline.scattering_angle_in_yz_plane( - incident_beam=incident_beam, - scattered_beam=scattered_beam, - wavelength=wavelength, - gravity=gravity, - ) - - sc.testing.assert_allclose(res, sc.scalar(polar, unit='rad')) - - @pytest.mark.parametrize('polar', [np.pi / 3, np.pi / 2, 2 * np.pi / 3, np.pi]) @pytest.mark.parametrize('x', [0.5, 11.4, -9.7]) def test_scattering_angle_in_yz_plane_does_not_depend_on_x(polar: float, x: float):