Gvec xy test

tests/test_transforms.py

@@ -1,3 +1,4 @@
+from collections import namedtuple
 import copy
 import json
 
@@ -8,6 +9,13 @@
 from common import convert_axis_angle_to_rmat
 
 
+_flds = [
+    "tvec_c", "tvec_d", "tvec_s", "rmat_c", "rmat_d", "rmat_s",
+    "gvec_c", "beam_vec", "xy"
+]
+GvecXYData = namedtuple("GvecData", _flds)
+
+
 def test_gvec_to_xy(test_data_dir):
     with open(test_data_dir / 'gvec_to_xy.json') as rf:
         test_data = json.load(rf)
@@ -36,3 +44,329 @@ def test_gvec_to_xy(test_data_dir):
 
         result = gvec_to_xy(**rotated_kwargs)
         assert np.allclose(result, expected)
+
+
+class TestGvecXY:
+    """Test gvec_to_xy and xy_to_gvec"""
+
+    # Base Case: sample and crystal align with lab frame;
+    # detector distance = 10; gvec = beam, no diffraction.
+    base = GvecXYData(
+        tvec_c=np.zeros(3),
+        tvec_d=np.array([0., 0., -10]),
+        tvec_s=np.zeros(3),
+        rmat_c=np.identity(3),
+        rmat_d=np.identity(3),
+        rmat_s=np.identity(3),
+        gvec_c=np.array([0, 0, -1.0]),
+        beam_vec=np.array([0, 0, -1.0]),
+        xy = np.array((np.nan, np.nan))
+    )
+
+    @staticmethod
+    def to_array(a):
+        return np.array(a, dtype=float)
+
+    @classmethod
+    def run_test(cls, prob):
+        """Run a test problem"""
+        xy_d = gvec_to_xy(
+            cls.to_array(prob.gvec_c),
+            cls.to_array(prob.rmat_d),
+            cls.to_array(prob.rmat_s),
+            cls.to_array(prob.rmat_c),
+            cls.to_array(prob.tvec_d),
+            cls.to_array(prob.tvec_s),
+            cls.to_array(prob.tvec_c),
+            beam_vec=cls.to_array(prob.beam_vec),
+        )
+        assert np.allclose(xy_d, prob.xy, equal_nan=True)
+
+    @staticmethod
+    def gvec_dvec(theta_deg, eta_deg):
+        """Return diffraction vector
+
+        PARAMETERS
+        ----------
+        theta_deg: float
+           diffraction angle theta, in degrees
+        eta_degrees: float
+           azimuthal angle, in degrees
+
+        RETURNS
+        -------
+        gvec: array(3):
+           diffraction vector relative to beam in -z-direction
+        dvec: array(3):
+           direction (unit vector) of diffracted beam
+        """
+        eta = np.radians(eta_deg)
+        ce, se = np.cos(eta), np.sin(eta)
+
+        # Diffraction vector makes angle theta with plane normal.
+        normal_ang = np.radians(90 - theta_deg)
+        cna, sna = np.cos(normal_ang), np.sin(normal_ang)
+        gvec = np.array((ce * sna, se * sna, cna))
+
+        # Diffracted beam direction makes angle 2-theta with beam.
+        two_theta = np.radians(2 * theta_deg)
+        c2t, s2t = np.cos(two_theta), np.sin(two_theta)
+        dvec = np.array((ce * s2t, se * s2t, -c2t))
+
+        return gvec, dvec
+
+    @staticmethod
+    def make_rmat(angle_deg, axis):
+        """Make a rotation matrix
+
+        PARAMETERS
+        ----------
+        angle: float
+           rotation angle in degrees
+        axis: array(3)
+           axis of rotation (not normalized to unit vector)
+
+        RETURNS
+        -------
+        array(3, 3):
+           rotation matrix
+        """
+        return convert_axis_angle_to_rmat(
+            np.array(axis), np.radians(angle_deg)
+        )
+
+    @staticmethod
+    def line_plane_intersect(p0, dvec, d0, nvec):
+        """Solve line/plane intersection
+
+        PARAMETERS
+        ----------
+        p0: array/3-tuple of floats
+            a point on the line
+        dvec: array/3-tuple of floats
+            direction (unit) vector of line
+        d0: array/3-tuple of floats
+            origin of detector
+        nvec: array/3-tuple of floats
+            normal vector to detector
+
+        RETURNS
+        -------
+        xyz: 3-tuple
+            intersection point
+
+        NOTES
+        -----
+        From geometry, t = (d0 - p0).n/(d.n), and x = p0 + td
+        """
+        p0_, d = np.array([p0, dvec])
+        t = np.dot((d0 - p0_), nvec)/np.dot(d, nvec)
+        x = p0_ + t * d
+        return x
+
+    @classmethod
+    def test_theta_eta(cls):
+        """Vary diffraction angle in simple case
+
+        TEST PARAMETERS
+        ---------------
+        theta_deg: float
+           diffraction angle in degrees
+        eta_deg: float, default=0
+           azimuthal angle in degrees
+        """
+        ThetaEta = namedtuple("ThetaEta", ["theta_deg", "eta_deg"])
+        nan_tests = [
+            ThetaEta(0, 0), ThetaEta(46, 0),
+        ]
+
+        tests = [
+            ThetaEta(10, 0), ThetaEta(44.9, 0),
+            ThetaEta(10, 45), ThetaEta(10, -45),
+        ]
+
+        p0_l = (0, 0, 0)
+        d0_l = cls.base.tvec_d
+        nv_l = (0, 0, 1)
+        for t in tests:
+            print("test: ", t)
+            gvec, dvec = cls.gvec_dvec(t.theta_deg, t.eta_deg)
+            det_x = cls.line_plane_intersect(
+                p0_l, dvec, d0_l, nv_l
+            )
+            if t.theta_deg <= 0:
+                answer = np.array((np.nan, np.nan))
+            else:
+                answer = det_x[:2]
+
+            cls.run_test(cls.base._replace(
+                gvec_c=gvec, xy=answer
+            ))
+
+    @classmethod
+    def test_beam(cls):
+        """Vary beam direction
+
+        The beam is rotated by a specified rotation matrix, determined
+        from an angle/axis pair.
+
+        TEST PARAMETERS
+        ---------------
+        aa: 2-tuple
+           angle/axis pairs
+        """
+        TestData = namedtuple("TestData", ["angle", "axis"])
+        tests = [
+            TestData(90, (1, 0, 0)),
+            TestData(90, (0, 1, 0)),
+            TestData(45, (1, 0, 0)),
+            TestData(45, (0, 0, 1))
+        ]
+
+        theta_deg, eta_deg = 5, 0
+        gvec, dvec = cls.gvec_dvec(theta_deg, eta_deg)
+
+        p0_l = (0, 0, 0)
+        d0_l = cls.base.tvec_d
+        nv_d = (0, 0, 1)
+        for t in tests:
+            print("test: ", t)
+            rmat_d = cls.make_rmat(t.angle, t.axis)
+            gvec_l = rmat_d @ gvec
+            dvec_l = rmat_d @ dvec
+            tvec_l = rmat_d @ cls.base.tvec_d
+            beam = rmat_d @ cls.base.beam_vec
+            det_x = cls.line_plane_intersect(
+                p0_l, dvec_l, tvec_l, rmat_d @ nv_d
+            )
+            x_d = rmat_d.T @ (det_x - tvec_l)
+
+            cls.run_test(cls.base._replace(
+                beam_vec=beam, tvec_d=tvec_l, gvec_c=gvec_l,
+                rmat_d=rmat_d, xy=x_d[:2],
+            ))
+
+    @classmethod
+    def test_translations(cls):
+        """Vary translation vectors
+
+        TEST PARAMETERS
+        ---------------
+        td: array(3)/3-tuple
+           detector translation in lab frame
+        ts:
+           sample translation in lab frame
+        tc:
+           crystal translation in sample frame
+        """
+        _flds = ["td", "ts", "tc"]
+        TvecData = namedtuple("TvecData", _flds)
+
+        tests = [
+            TvecData(td=(0, 0, -10), ts=(0, 0, 0), tc=(0, 0, 0)),
+            TvecData(td=(0, 0, -10), ts=(0, 0, 7), tc=(0, 0, 0)),
+            TvecData(td=(2, 3, -10), ts=(0, 0, 0), tc=(0, 0, -3)),
+            TvecData(td=(0, 0, -10), ts=(0, 1, 2), tc=(2, 3, -4)),
+        ]
+
+        nv_l = (0, 0, 1)
+        theta_deg, eta_deg = 7.5, 0
+        gvec, dvec = cls.gvec_dvec(theta_deg, eta_deg)
+        for t in tests:
+            print("test data:\n", t)
+            td = np.array(t.td)
+            ts = np.array(t.ts)
+            tc = np.array(t.tc)
+            p0_l = ts + tc
+            det_x = cls.line_plane_intersect(
+                p0_l, dvec, td, nv_l
+            )
+            xy = det_x[:2] - td[:2]
+
+            cls.run_test(cls.base._replace(
+                gvec_c=gvec, tvec_d=td, tvec_s=ts, tvec_c=tc,
+                xy=xy
+            ))
+
+    @classmethod
+    def test_detector_orientation(cls):
+        """Vary detector orientation
+
+        TEST PARAMETERS
+        ---------------
+        rd: 2-tuple
+           angle/axis pairs
+        """
+        RmatData = namedtuple("RmatData", ["angle_deg", "axis"])
+        nan_tests = [
+            RmatData(90, (0, 1, 0)),
+        ]
+        tests = [
+            RmatData(90, (1, 0, 0)),
+            RmatData(45, (1, 0, 0)),
+            RmatData(45, (0, 0, 1))
+        ]
+
+        theta_deg, eta_deg = 5, 0
+        gvec, dvec = cls.gvec_dvec(theta_deg, eta_deg)
+
+        p0_l = (0, 0, 0)
+        d0_l = cls.base.tvec_d
+        nv_l = (0, 0, 1)
+        for t in tests:
+            print(t)
+            tvec_l = -cls.base.tvec_d
+            rmat_d = cls.make_rmat(t.angle_deg, t.axis)
+            det_x = cls.line_plane_intersect(
+                p0_l, dvec, d0_l, rmat_d @ nv_l
+            )
+            x_d = rmat_d.T @ (det_x - d0_l)
+            cls.run_test(cls.base._replace(
+                gvec_c=gvec, rmat_d=rmat_d, xy=x_d[:2]
+            ))
+
+    @classmethod
+    def test_sample_crystal_orientations(cls):
+        """Vary sample and crystal orientations
+
+        TEST PARAMETERS
+        ----------
+        data: named_tuple
+           angle/axis pairs for sample and crystal rotation
+        """
+        _flds = ["aa_s", "aa_c"]
+        SampleCrystalData = namedtuple("SampleCrystalData", _flds)
+
+        ex, ey, ez = (1, 0, 0), (0, 1, 0), (0, 0, 1)
+        ang = 10
+        tests = [
+            SampleCrystalData((0, ex), (0, ex)),
+            SampleCrystalData((ang, ex), (-ang, ex)),
+            SampleCrystalData((ang, ey), (-ang, ey)),
+            SampleCrystalData((ang, ez), (-ang, ez)),
+            SampleCrystalData((ang, ex), (0, ex)),
+            SampleCrystalData((ang, ey), (0, ey)),
+            SampleCrystalData((ang, ez), (0, ez)),
+            SampleCrystalData((0, ex), (ang, ex)),
+            SampleCrystalData((0, ey), (ang, ey)),
+            SampleCrystalData((0, ez), (ang, ez)),
+        ]
+
+        theta_deg, eta_deg = 5, 90
+        gvec_l, dvec_l = cls.gvec_dvec(theta_deg, eta_deg)
+        for t in tests:
+            print(t)
+            ang_s, ax_s = t.aa_s
+            rmat_s = cls.make_rmat(ang_s, ax_s)
+            ang_c, ax_c = t.aa_c
+            rmat_c = cls.make_rmat(ang_c, ax_c)
+
+            gvec_c = rmat_c.T @ rmat_s.T @ gvec_l
+            det_x = cls.line_plane_intersect(
+                (0, 0, 0), dvec_l, cls.base.tvec_d, (0, 0, 1)
+            )
+            xy = det_x[:2]
+
+            cls.run_test(cls.base._replace(
+                gvec_c=gvec_c, rmat_s=rmat_s, rmat_c=rmat_c, xy=xy
+            ))