Merge pull request #633 from scipp/fix-coord-conv-position-unit-mismatch-exception

Fix coord conversion position unit mismatch exceptions

Author: SimonHeybrock

src/scippneutron/conversion/beamline.py

@@ -117,6 +117,11 @@
 from .._utils import elem_dtype, elem_unit
 
 
+def _canonical_length(var: VariableLike) -> VariableLike:
+    """Convert a beam length to the canonical unit (meter)."""
+    return var.to(unit='m', copy=False)
+
+
 def L1(*, incident_beam: VariableLike) -> VariableLike:
     """Compute the length of the incident beam.
 
@@ -141,7 +146,7 @@ def L1(*, incident_beam: VariableLike) -> VariableLike:
     straight_incident_beam:
         Compute the incident beam for a straight beamline.
     """
-    return sc.norm(incident_beam)
+    return _canonical_length(sc.norm(incident_beam))
 
 
 def L2(*, scattered_beam: VariableLike) -> VariableLike:
@@ -168,7 +173,7 @@ def L2(*, scattered_beam: VariableLike) -> VariableLike:
     straight_scattered_beam:
         Compute the scattered beam for a straight beamline.
     """
-    return sc.norm(scattered_beam)
+    return _canonical_length(sc.norm(scattered_beam))
 
 
 def straight_incident_beam(
@@ -196,7 +201,7 @@ def straight_incident_beam(
     :
         ``incident_beam``
     """
-    return sample_position - source_position
+    return _canonical_length(sample_position) - _canonical_length(source_position)
 
 
 def straight_scattered_beam(
@@ -223,7 +228,7 @@ def straight_scattered_beam(
     :
         ``scattered_beam``
     """
-    return position - sample_position
+    return _canonical_length(position) - _canonical_length(sample_position)
 
 
 def total_beam_length(*, L1: VariableLike, L2: VariableLike) -> VariableLike:
@@ -247,7 +252,7 @@ def total_beam_length(*, L1: VariableLike, L2: VariableLike) -> VariableLike:
     :
         :math:`L_\\mathsf{total}`
     """
-    return L1 + L2
+    return _canonical_length(L1) + _canonical_length(L2)
 
 
 def total_straight_beam_length_no_scatter(
@@ -274,7 +279,7 @@ def total_straight_beam_length_no_scatter(
     :
         :math:`L_\\mathsf{total}`
     """
-    return sc.norm(position - source_position.to(unit=position.unit, copy=False))
+    return sc.norm(_canonical_length(position) - _canonical_length(source_position))
 
 
 def two_theta(
@@ -324,6 +329,8 @@ def two_theta(
     # and 'it never errs by more than a modest multiple of epsilon'
     # where 'epsilon is the roundoff threshold for individual arithmetic
     # operations'.
+    incident_beam = _canonical_length(incident_beam)
+    scattered_beam = _canonical_length(scattered_beam)
     b1 = incident_beam / L1(incident_beam=incident_beam)
     b2 = scattered_beam / L2(scattered_beam=scattered_beam)
 
@@ -543,6 +550,8 @@ def scattering_angles_with_gravity(
         Ignores the ``x`` component when computing ``theta``.
         This is used in reflectometry.
     """
+    incident_beam = _canonical_length(incident_beam)
+    scattered_beam = _canonical_length(scattered_beam)
     if sc.any(
         abs(sc.dot(gravity, incident_beam))
         > sc.scalar(1e-10, unit=incident_beam.unit) * sc.norm(gravity)

tests/conversion/beamline_conversions_test.py

@@ -35,7 +35,9 @@ def test_straight_incident_beam():
     )
     sc.testing.assert_allclose(
         incident_beam,
-        sc.vectors(dims=['siti'], values=[[-33, -75, -17], [-23, -5, -47]], unit='mm'),
+        sc.vectors(
+            dims=['siti'], values=[[-33, -75, -17], [-23, -5, -47]], unit='mm'
+        ).to(unit='m'),
     )
 
 
@@ -47,7 +49,9 @@ def test_straight_scattered_beam():
     )
     sc.testing.assert_allclose(
         scattered_beam,
-        sc.vectors(dims=['on'], values=[[5.6, 2.4, 8.8], [7.6, 2.4, 0.8]], unit='km'),
+        sc.vectors(
+            dims=['on'], values=[[5.6, 2.4, 8.8], [7.6, 2.4, 0.8]], unit='km'
+        ).to(unit='m'),
     )
 
 
@@ -57,7 +61,10 @@ def test_L1():
     )
     L1 = beamline.L1(incident_beam=incident_beam)
     sc.testing.assert_allclose(
-        L1, sc.array(dims=['inc'], values=[1.870828693386, 1.122497216032], unit='um')
+        L1,
+        sc.array(dims=['inc'], values=[1.870828693386, 1.122497216032], unit='um').to(
+            unit='m'
+        ),
     )
 
 
@@ -68,7 +75,9 @@ def test_L2():
     L2 = beamline.L2(scattered_beam=scattered_beam)
     sc.testing.assert_allclose(
         L2,
-        sc.array(dims=['scat'], values=[41.158231254513, 146.321563687653], unit='am'),
+        sc.array(
+            dims=['scat'], values=[41.158231254513, 146.321563687653], unit='am'
+        ).to(unit='m'),
     )
 
 
@@ -78,7 +87,9 @@ def test_total_beam_length():
     Ltotal = beamline.total_beam_length(L1=L1, L2=L2)
     sc.testing.assert_allclose(
         Ltotal,
-        sc.array(dims=['secondary'], values=[8.27559, 47.134, 1004.134], unit='cm'),
+        sc.array(dims=['secondary'], values=[8.27559, 47.134, 1004.134], unit='cm').to(
+            unit='m'
+        ),
     )
 
 
@@ -672,7 +683,7 @@ def test_scattering_angles_with_gravity_supports_mismatching_units(g: list[float
     sc.testing.assert_allclose(
         res['two_theta'], expected['two_theta'], rtol=sc.scalar(1e-6)
     )
-    sc.testing.assert_allclose(res['phi'], expected['phi'])
+    sc.testing.assert_allclose(res['phi'], expected['phi'], rtol=sc.scalar(1e-6))
 
 
 @given(rotation=rotations())