add functions to allow conversions between d and tth, and d and q

news/xtype.rst

@@ -0,0 +1,23 @@
+**Added:**
+
+* functions to allow conversions between d and tth, and d and q.
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

src/diffpy/utils/scattering_objects/diffraction_objects.py

@@ -10,6 +10,8 @@
 DQUANTITIES = ["d", "dspace"]
 XQUANTITIES = ANGLEQUANTITIES + DQUANTITIES + QQUANTITIES
 XUNITS = ["degrees", "radians", "rad", "deg", "inv_angs", "inv_nm", "nm-1", "A-1"]
+QMAX = 40
+DMAX = 100
 
 x_grid_emsg = (
     "objects are not on the same x-grid. You may add them using the self.add method "
@@ -283,7 +285,7 @@ def insert_scattering_quantity(
         elif xtype.lower() in ANGLEQUANTITIES:
             self.on_tth = [np.array(xarray), np.array(yarray)]
         elif xtype.lower() in DQUANTITIES:
-            self.on_tth = [np.array(xarray), np.array(yarray)]
+            self.on_d = [np.array(xarray), np.array(yarray)]
         self.set_all_arrays()
 
     def q_to_tth(self):
@@ -305,7 +307,7 @@ def q_to_tth(self):
         Parameters
         ----------
         q : array
-            An array of :math:`q` values
+            The array of :math:`q` values
 
         wavelength : float
             Wavelength of the incoming x-rays
@@ -315,13 +317,15 @@ def q_to_tth(self):
         Returns
         -------
         two_theta : array
-            An array of :math:`2\theta` values in radians
+            The array of :math:`2\theta` values in radians
         """
         q = self.on_q[0]
         q = np.asarray(q)
         wavelength = float(self.wavelength)
         pre_factor = wavelength / (4 * np.pi)
-        return np.rad2deg(2.0 * np.arcsin(q * pre_factor))
+        # limit q * pre_factor to 1 to avoid invalid input to arcsin
+        arcsin_value = np.where(q * pre_factor <= 1, q * pre_factor, 1)
+        return np.rad2deg(2.0 * np.arcsin(arcsin_value))
 
     def tth_to_q(self):
         r"""
@@ -344,7 +348,7 @@ def tth_to_q(self):
         Parameters
         ----------
         two_theta : array
-            An array of :math:`2\theta` values in units of degrees
+            The array of :math:`2\theta` values in units of degrees
 
         wavelength : float
             Wavelength of the incoming x-rays
@@ -354,26 +358,126 @@ def tth_to_q(self):
         Returns
         -------
         q : array
-            An array of :math:`q` values in the inverse of the units
+            The array of :math:`q` values in the inverse of the units
             of ``wavelength``
         """
         two_theta = np.asarray(np.deg2rad(self.on_tth[0]))
         wavelength = float(self.wavelength)
         pre_factor = (4 * np.pi) / wavelength
         return pre_factor * np.sin(two_theta / 2)
 
+    def q_to_d(self):
+        r"""
+        Helper function to convert q to d using :math:`d = \frac{2 \pi}{q}`, set dmax to DMAX=100
+
+        Parameters
+        ----------
+        q : array
+            The array of :math:`q` values
+
+        Returns
+        -------
+        d : array
+            The array of :math:`d` values in the inverse of the units of ``wavelength``
+        """
+        q = np.asarray(self.on_q[0])
+        d = np.where(q != 0, (2 * np.pi) / q, np.inf)
+        d = np.minimum(d, DMAX)
+        return d[::-1]
+
+    def d_to_q(self):
+        r"""
+        Helper function to convert d to q using :math:`q = \frac{2 \pi}{d}`, set qmax to QMAX=40
+
+        Parameters
+        ----------
+        d : array
+            The array of :math:`d` values
+
+        Returns
+        -------
+        q : array
+            The array of :math:`q` values in the inverse of the units of ``wavelength``
+        """
+        d = np.asarray(self.on_d[0])
+        q = np.where(d != 0, (2 * np.pi) / d, np.inf)
+        q = np.minimum(q, QMAX)
+        return q[::-1]
+
+    def tth_to_d(self):
+        r"""
+        Helper function to convert two-theta to d
+
+        uses the formula .. math:: d = \frac{\lambda}{2 \sin\left(\frac{2\theta}{2}\right)}, set dmax to DMAX=100
+
+        Parameters
+        ----------
+        two_theta : array
+            The array of :math:`2\theta` values in units of degrees
+
+        wavelength : float
+            Wavelength of the incoming x-rays
+
+        Returns
+        -------
+        d : array
+            The array of :math:`d` values in the inverse of the units
+            of ``wavelength``
+        """
+        two_theta = np.asarray(np.deg2rad(self.on_tth[0]))
+        wavelength = float(self.wavelength)
+        sin_two_theta = np.sin(two_theta / 2)
+        d = np.where(sin_two_theta != 0, wavelength / (2 * sin_two_theta), np.inf)
+        d = np.minimum(d, DMAX)
+        return d[::-1]
+
+    def d_to_tth(self):
+        r"""
+        Helper function to convert d to two-theta
+
+        uses the formula .. math:: 2\theta = 2 \arcsin\left(\frac{\lambda}{2d}\right), set tth to 180 when d=0
+
+        Parameters
+        ----------
+        d : array
+            The array of :math:`d` values
+
+        wavelength : float
+            Wavelength of the incoming x-rays
+
+        Returns
+        -------
+        two_theta : array
+            The array of :math:`2\theta` values in radians
+        """
+        d = np.asarray(self.on_d[0])
+        wavelength = float(self.wavelength)
+        pre_factor = np.where(d != 0, wavelength / (2 * d), 1)
+        return np.rad2deg(2.0 * np.arcsin(pre_factor))[::-1]
+
     def set_all_arrays(self):
         master_array, xtype = self._get_original_array()
         if xtype == "q":
             self.on_tth[0] = self.q_to_tth()
             self.on_tth[1] = master_array[1]
+            self.on_d[0] = self.q_to_d()
+            self.on_d[1] = master_array[1]
         if xtype == "tth":
             self.on_q[0] = self.tth_to_q()
             self.on_q[1] = master_array[1]
+            self.on_d[0] = self.tth_to_d()
+            self.on_d[1] = master_array[1]
+        if xtype == "d":
+            self.on_tth[0] = self.d_to_tth()
+            self.on_tth[1] = master_array[1]
+            self.on_q[0] = self.d_to_q()
+            self.on_q[1] = master_array[1]
         self.tthmin = self.on_tth[0][0]
         self.tthmax = self.on_tth[0][-1]
         self.qmin = self.on_q[0][0]
         self.qmax = self.on_q[0][-1]
+        self.dmin = self.on_d[0][0]
+        self.dmax = self.on_d[0][-1]
 
     def _get_original_array(self):
         if self.input_xtype in QQUANTITIES:

tests/test_diffraction_objects.py

@@ -231,6 +231,95 @@ def test_diffraction_objects_equality(inputs1, inputs2, expected):
     assert (diffraction_object1 == diffraction_object2) == expected
 
 
+def test_q_to_tth():
+    actual = Diffraction_object(wavelength=4 * np.pi)
+    setattr(actual, "on_q", [[0, 0.2, 0.4, 0.6, 0.8, 1, 40, 60], [1, 2, 3, 4, 5, 6, 7, 8]])
+    actual_tth = actual.q_to_tth()
+    # expected tth values are 2 * arcsin(q)
+    # when q gets large, we set tth = 180
+    # we allow q to exceed QMAX for user inputs
+    expected_tth = [0, 23.07392, 47.15636, 73.73980, 106.26020, 180, 180, 180]
+    assert np.allclose(actual_tth, expected_tth)
+
+
+def test_tth_to_q():
+    actual = Diffraction_object(wavelength=4 * np.pi)
+    setattr(actual, "on_tth", [[0, 30, 60, 90, 120, 180], [1, 2, 3, 4, 5, 6]])
+    actual_q = actual.tth_to_q()
+    # expected q vales are sin15, sin30, sin45, sin60, sin90
+    expected_q = [0, 0.258819, 0.5, 0.707107, 0.866025, 1]
+    assert np.allclose(actual_q, expected_q)
+
+
+def test_q_to_d():
+    actual = Diffraction_object()
+    setattr(actual, "on_q", [[0, np.pi, 2 * np.pi, 3 * np.pi, 4 * np.pi, 5 * np.pi], [1, 2, 3, 4, 5, 6]])
+    actual_d = actual.q_to_d()
+    # expected d values are DMAX=100, 2/1, 2/2, 2/3, 2/4, 2/5, in reverse order
+    expected_d = [0.4, 0.5, 0.66667, 1, 2, 100]
+    assert np.allclose(actual_d, expected_d)
+
+
+def test_d_to_q():
+    actual = Diffraction_object()
+    setattr(actual, "on_d", [[0, np.pi, 2 * np.pi, 3 * np.pi, 4 * np.pi, 5 * np.pi], [1, 2, 3, 4, 5, 6]])
+    actual_q = actual.d_to_q()
+    # expected q values are QMAX=40, 2/1, 2/2, 2/3, 2/4, 2/5, in reverse order
+    expected_q = [0.4, 0.5, 0.66667, 1, 2, 40]
+    assert np.allclose(actual_q, expected_q)
+
+
+def test_tth_to_d():
+    actual = Diffraction_object(wavelength=2)
+    setattr(actual, "on_tth", [[0, 30, 60, 90, 120, 180], [1, 2, 3, 4, 5, 6]])
+    actual_d = actual.tth_to_d()
+    # expected d values are DMAX=100, 1/sin15, 1/sin30, 1/sin45, 1/sin60, 1/sin90, in reverse order
+    expected_d = [1, 1.1547, 1.41421, 2, 3.8637, 100]
+    assert np.allclose(actual_d, expected_d)
+
+
+def test_d_to_tth():
+    actual = Diffraction_object(wavelength=2)
+    setattr(actual, "on_d", [[0, 2, 4, 6, 8, 100, 200], [1, 2, 3, 4, 5, 6, 7]])
+    actual_tth = actual.d_to_tth()
+    # expected tth values are 2*arcsin(1/d), in reverse order
+    # when d is 0, we set tth to 180
+    # we allow d to exceed DMAX=100 for user inputs
+    expected_tth = [0.57296, 1.14593, 14.36151, 19.18814, 28.95502, 60, 180]
+    assert np.allclose(actual_tth, expected_tth)
+
+
+params_array = [
+    (["tth", "on_tth", [30, 60, 90, 120, 150], [1, 2, 3, 4, 5]]),
+    (["q", "on_q", [1.626208, 3.141593, 4.442883, 5.441398, 6.069091], [1, 2, 3, 4, 5]]),
+    (["d", "on_d", [1.035276, 1.154701, 1.414214, 2, 3.863703], [1, 2, 3, 4, 5]]),
+]
+
+
+@pytest.mark.parametrize("inputs", params_array)
+def test_set_all_arrays(inputs):
+    input_xtype, on_xtype, xarray, yarray = inputs
+    expected_values = {
+        "on_tth": [np.array([30, 60, 90, 120, 150]), np.array([1, 2, 3, 4, 5])],
+        "on_q": [np.array([1.626208, 3.141593, 4.442883, 5.441398, 6.069091]), np.array([1, 2, 3, 4, 5])],
+        "on_d": [np.array([1.035276, 1.154701, 1.414214, 2, 3.863703]), np.array([1, 2, 3, 4, 5])],
+        "tthmin": 30,
+        "tthmax": 150,
+        "qmin": 1.626208,
+        "qmax": 6.069091,
+        "dmin": 1.035276,
+        "dmax": 3.863703,
+    }
+
+    actual = Diffraction_object(wavelength=2)
+    setattr(actual, "input_xtype", input_xtype)
+    setattr(actual, on_xtype, [xarray, yarray])
+    actual.set_all_arrays()
+    for attr, expected in expected_values.items():
+        actual_value = getattr(actual, attr)
+        assert np.allclose(actual_value, expected)
+
+
 def test_dump(tmp_path, mocker):
     x, y = np.linspace(0, 5, 6), np.linspace(0, 5, 6)
     directory = Path(tmp_path)