transform_coords building blocks for gravity correction in theta vs wavelength transformations

src/ess/reflectometry/transform_coords.py

@@ -0,0 +1,67 @@
+import scipp as sc
+import numpy as np
+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)
+
+
+def to_y_dash(wavelength, scattered_beam, vertical_unit_vector,
+              forward_beam_unit_vector):
+    velocity_sq = to_velocity(wavelength)
+    velocity_sq *= velocity_sq
+    g_v = sc.norm(vertical_unit_vector * g)
+    # dy due to gravity = -0.5gt^2 = -0.5g(dz/dv)^2
+    # therefore y'(z) = dy/dz - 0.5g.dz/dv^2 / dz
+    forward = sc.dot(scattered_beam, forward_beam_unit_vector)
+    vertical = sc.dot(scattered_beam, vertical_unit_vector)
+    return (-0.5 * g_v * forward / velocity_sq) + (vertical / forward)
+
+
+def incident_beam(sample_position, source_position):
+    return sample_position - source_position
+
+
+def scattered_beam(position, sample_position):
+    return position - sample_position
+
+
+def L1(incident_beam):
+    return sc.norm(incident_beam)
+
+
+def L2(scattered_beam):
+    return sc.norm(scattered_beam)
+
+
+def _angle(a, b):
+    return sc.acos(sc.dot(a, b) / (sc.norm(a) * sc.norm(b)))
+
+
+def to_scattering_angle(w_norm, wavelength, detector_id, sample_position, incident_beam,
+                        scattered_beam):
+    w_norm = w_norm / sc.norm(w_norm)
+    incident_beam_norm = incident_beam / sc.norm(incident_beam)
+    scattered_beam_norm = scattered_beam / sc.norm(scattered_beam)
+    # vector pointing along surface in forward beam direction
+    surface = sc.cross(w_norm, sc.cross(incident_beam_norm, scattered_beam_norm))
+    # Assume specular reflection. And reflect incident beam through surface
+    reflection = incident_beam - 2.0 * sc.dot(incident_beam, surface) * surface
+    forward_beam_direction = incident_beam - reflection
+    # For a specular reflection, this would be basis aligned
+    forward_beam_direction /= sc.norm(forward_beam_direction)
+    # Account for non-specular scattering
+    forward_beam_direction = sc.vector(value=np.round(forward_beam_direction.value),
+                                       unit=forward_beam_direction.unit)
+    # Vertical direction
+    vertical_direction = sc.vector(value=np.round(w_norm.value), unit=w_norm.unit)
+    y_dash = to_y_dash(wavelength, scattered_beam, vertical_direction,
+                       forward_beam_direction)
+    start = sc.dot(vertical_direction, sample_position)
+    height = y_dash * sc.dot(forward_beam_direction, scattered_beam) + start
+
+    w = _angle(w_norm, surface) - sc.scalar(value=np.pi / 2, unit=sc.units.rad)
+    return sc.atan2(y=height, x=sc.dot(forward_beam_direction, incident_beam)) - w

tests/reflectometry/transform_coords_test.py

@@ -0,0 +1,138 @@
+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)
+    scattered_beam = detector_position - sample_position
+    vertical_unit_vector = sc.vector(value=[0, 1, 0])
+    forward_beam_unit_vector = sc.vector(value=[0, 0, 1])
+
+    # 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,
+                                      scattered_beam=scattered_beam,
+                                      vertical_unit_vector=vertical_unit_vector,
+                                      forward_beam_unit_vector=forward_beam_unit_vector)
+
+    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)
+    scattered_beam = detector_position - sample_position
+    vertical_unit_vector = sc.vector(value=[0, 1, 0])
+    fwd_beam_unit_vector = sc.vector(value=[0, 0, 1])
+
+    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
+    transform_args = {
+        "wavelength": wav,
+        "scattered_beam" : scattered_beam,
+        "vertical_unit_vector" : vertical_unit_vector,
+        "forward_beam_unit_vector" : fwd_beam_unit_vector
+    }
+    grad = transform_coords.to_y_dash(**transform_args)
+    assert sc.less(grad, 1 * sc.units.one).value
+
+    vel *= 2
+    transform_args["wavelength"]\
+        = sc.to_unit(h / (vel * neutron_mass), unit=sc.units.angstrom)
+    grad_fast = transform_coords.to_y_dash(**transform_args)
+    # 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():
+    source_position = sc.vector(value=[0, 1, -1], unit=sc.units.m)
+    sample_position = sc.vector(value=[0, 0, 0], unit=sc.units.m)
+    detector_position = sc.vector(value=[0, 1, 1], unit=sc.units.m)
+    incident_beam = sample_position - source_position
+    scattered_beam = detector_position - sample_position
+    no_gravity_angle = _angle(scattered_beam, incident_beam)
+
+    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_norm=sc.vector(value=[0, 1, 0]),
+                                                 wavelength=wav,
+                                                 detector_id=None,
+                                                 sample_position=sample_position,
+                                                 incident_beam=incident_beam,
+                                                 scattered_beam=scattered_beam)
+    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]), incident_beam) / 2.0
+    assert sc.isclose(angle, expected).value
+
+
+def test_scattering_angle_xzy():
+    # Same as previous but we define forward beam direction to be +x
+    # up direction to be z (gravity therefore acts in -z)
+    # perpendicular direction to be y, as in w is rotation around y
+
+    source_position = sc.vector(value=[-1, 0, 1], unit=sc.units.m)
+    sample_position = sc.vector(value=[0, 0, 0], unit=sc.units.m)
+    detector_position = sc.vector(value=[1, 0, 1], unit=sc.units.m)
+    incident_beam = sample_position - source_position
+    scattered_beam = detector_position - sample_position
+
+    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_norm=sc.vector(value=[0, 0, 1]),
+                                                 wavelength=wav,
+                                                 detector_id=None,
+                                                 sample_position=sample_position,
+                                                 incident_beam=incident_beam,
+                                                 scattered_beam=scattered_beam)
+
+    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]), incident_beam) / 2.0
+    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, 1, -1], 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_norm"] = sc.vector(value=[0, 1, 0], unit=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,
+                "incident_beam": transform_coords.incident_beam,
+                "scattered_beam": transform_coords.scattered_beam}
+    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