transform_coords building blocks for gravity correction in theta vs wavelength transformations

src/ess/reflectometry/transform_coords.py

@@ -0,0 +1,28 @@
+import scipp as sc
+from scipp.constants import neutron_mass, h, g
+
+
+def to_velocity(wavelength):
+    return sc.to_unit(h / (wavelength * neutron_mass), sc.units.m / sc.units.s,
+                      copy=False)
+
+
+# Derivative of y with respect to z
+def to_y_dash(wavelength, sample_position, detector_position):
+    diff = (detector_position - sample_position)
+    velocity_sq = to_velocity(wavelength)
+    velocity_sq *= velocity_sq
+    gy = sc.vector(value=[0, -1, 0]) * g
+    # dy due to gravity = -0.5gt^2 = -0.5g(dz/dv)^2
+    # therefore y'(z) = dy/dz - 0.5g.dz/dv^2 / dz
+    return (-0.5 * sc.norm(gy)
+            * diff.fields.z / velocity_sq) + (diff.fields.y / diff.fields.z)
+
+
+def to_scattering_angle(w, wavelength, detector_id, position, sample_position):
+    z_origin = sample_position.fields.z
+    y_origin = sample_position.fields.y
+    z_measured = position.fields.z
+    y_dash = to_y_dash(wavelength, sample_position, position)
+    height = y_dash * (z_measured - z_origin) + y_origin
+    return sc.atan2(y=height, x=z_measured) - w

tests/reflectometry/transform_coords_test.py

@@ -0,0 +1,109 @@
+import scipp as sc
+from scipp.constants import g, h, neutron_mass
+from ess.reflectometry import transform_coords
+import numpy as np
+
+
+def test_y_dash_for_gravitational_effect():
+    sample_position = sc.vector(value=[0, 0, 0], unit=sc.units.m)
+    detector_position = sc.vector(value=[0, 0.5, 1], unit=sc.units.m)
+
+    # Approximate cold-neutron velocities
+    vel = 1000 * (sc.units.m / sc.units.s)
+    wav = sc.to_unit(h / (vel * neutron_mass),
+                     unit=sc.units.angstrom)
+
+    grad = transform_coords.to_y_dash(wavelength=wav, sample_position=sample_position,
+                                      detector_position=detector_position)
+
+    scattered_beam = detector_position - sample_position
+    no_gravity_grad = scattered_beam.fields.y / scattered_beam.fields.z
+    gravity_effect_grad = (-0.5 * g * scattered_beam.fields.z / (vel * vel))
+    assert sc.isclose(grad, no_gravity_grad + gravity_effect_grad).value
+
+
+def test_y_dash_with_different_velocities():
+    sample_position = sc.vector(value=[0, 0, 0], unit=sc.units.m)
+    detector_position = sc.vector(value=[0, 1, 1], unit=sc.units.m)
+
+    vel = 1000 * (sc.units.m / sc.units.s)
+    wav = sc.to_unit(h / (vel * neutron_mass),
+                     unit=sc.units.angstrom)
+
+    # In this setup the faster the neutrons the closer d'y(z) tends to 1.0
+    grad = transform_coords.to_y_dash(wavelength=wav,
+                                      sample_position=sample_position,
+                                      detector_position=detector_position)
+    assert sc.less(grad, 1 * sc.units.one).value
+
+    vel *= 2
+    wav = sc.to_unit(h / (vel * neutron_mass),
+                     unit=sc.units.angstrom)
+    grad_fast = transform_coords.to_y_dash(wavelength=wav,
+                                           sample_position=sample_position,
+                                           detector_position=detector_position)
+    # Testing that gravity has greater influence on slow neutrons.
+    assert sc.less(grad, grad_fast).value
+
+
+def _angle(a, b):
+    return sc.acos(sc.dot(a, b) / (sc.norm(a) * sc.norm(b)))
+
+
+def test_scattering_angle():
+    sample_position = sc.vector(value=[0, 0, 0], unit=sc.units.m)
+    detector_position = sc.vector(value=[0, 1, 1], unit=sc.units.m)
+    scattered_beam = detector_position - sample_position
+    beam_direction = sc.vector(value=[0, 0, 1], unit=sc.units.m)
+    no_gravity_angle = _angle(scattered_beam, beam_direction)
+
+    vel = 1000 * (sc.units.m / sc.units.s)
+    wav = sc.to_unit(h / (vel * neutron_mass),
+                     unit=sc.units.angstrom)
+
+    angle = transform_coords.to_scattering_angle(w=0 * sc.units.rad,
+                                                 wavelength=wav,
+                                                 detector_id=None,
+                                                 position=detector_position,
+                                                 sample_position=sample_position)
+    assert sc.less(angle, no_gravity_angle).value
+
+    gravity_shift_y = -0.5 * g * (scattered_beam.fields.z ** 2 / vel ** 2)
+    expected = _angle(scattered_beam + gravity_shift_y
+                      * sc.vector(value=[0, 1, 0]), beam_direction)
+    assert sc.isclose(angle, expected).value
+
+
+def test_det_wavelength_to_wavelength_scattering_angle():
+    # comparible with cold-neutrons from moderator
+    vel = 2000 * (sc.units.m / sc.units.s)
+    wav = sc.to_unit(h / (vel * neutron_mass),
+                     unit=sc.units.angstrom)
+    sample_position = sc.vector(value=[0, 0, 0], unit=sc.units.m)
+    source_position = sc.vector(value=[0, 0, -10], unit=sc.units.m)
+    detector_position = sc.vector(value=[0, 1, 1], unit=sc.units.m)
+
+    coords = {}
+    coords["sample_position"] = sample_position
+    coords["source_position"] = source_position
+    coords["position"] = detector_position
+    coords["wavelength"] = wav
+    coords["w"] = 0.0 * sc.units.rad
+    coords["detector_id"] = 0.0 * sc.units.one
+    measurement = sc.DataArray(data=1.0 * sc.units.one, coords=coords)
+
+    settings = {'scattering_angle': transform_coords.to_scattering_angle}
+    transformed = sc.transform_coords(x=measurement,
+                                      coords=['wavelength', 'scattering_angle'],
+                                      graph=settings)
+    assert sc.isclose(transformed.coords['scattering_angle'],
+                      (np.pi / 4) * sc.units.rad, atol=1e-4 * sc.units.rad).value
+
+    # We now check the sample angle. Setting to previous final scattering angle
+    # should yield a scattering angle of 0.
+    measurement.coords["w"] = transformed.coords['scattering_angle']
+    transformed = sc.transform_coords(x=measurement,
+                                      coords=['wavelength', 'scattering_angle'],
+                                      graph=settings)
+    assert sc.isclose(transformed.coords['scattering_angle'],
+                      0.0 * sc.units.rad).value