diff --git a/news/xtype.rst b/news/xtype.rst new file mode 100644 index 00000000..dbe976ec --- /dev/null +++ b/news/xtype.rst @@ -0,0 +1,23 @@ +**Added:** + +* functions to allow conversions between d and tth, and d and q. + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/src/diffpy/utils/scattering_objects/diffraction_objects.py b/src/diffpy/utils/scattering_objects/diffraction_objects.py index bd5b16db..6b6e490b 100644 --- a/src/diffpy/utils/scattering_objects/diffraction_objects.py +++ b/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,12 +317,16 @@ 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) + if np.any(np.abs(q * pre_factor) > 1): + raise ValueError( + "Invalid input for arcsin: some values exceed the range [-1, 1]. " "Check wavelength or q values." + ) return np.rad2deg(2.0 * np.arcsin(q * pre_factor)) def tth_to_q(self): @@ -344,7 +350,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 +360,134 @@ 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])) + if np.any(two_theta > np.pi): + raise ValueError("Two theta exceeds 180 degrees.") 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])) + if np.any(two_theta > np.pi): + raise ValueError("Two theta exceeds 180 degrees.") + 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 = wavelength / (2 * d) + if np.any(np.abs(pre_factor) > 1): + raise ValueError( + "Invalid input for arcsin: some values exceed the range [-1, 1]. " "Check wavelength or d values." + ) + 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: diff --git a/tests/test_diffraction_objects.py b/tests/test_diffraction_objects.py index a155b95e..dcaae022 100644 --- a/tests/test_diffraction_objects.py +++ b/tests/test_diffraction_objects.py @@ -231,6 +231,128 @@ def test_diffraction_objects_equality(inputs1, inputs2, expected): assert (diffraction_object1 == diffraction_object2) == expected +def test_q_to_tth(): + # valid q values including edge cases when q=0 or 1, expected tth values are 2 * arcsin(q) + actual = Diffraction_object(wavelength=4 * np.pi) + setattr(actual, "on_q", [[0, 0.2, 0.4, 0.6, 0.8, 1], [1, 2, 3, 4, 5, 6]]) + actual_tth = actual.q_to_tth() + expected_tth = [0, 23.07392, 47.15636, 73.73980, 106.26020, 180] + assert np.allclose(actual_tth, expected_tth) + + +def test_q_to_tth_bad(): + # invalid q values when arcsin value is not in the range of [-1, 1] + actual = Diffraction_object(wavelength=4 * np.pi) + setattr(actual, "on_q", [[0.6, 0.8, 1, 1.2], [1, 2, 3, 4]]) + with pytest.raises(ValueError): + actual.q_to_tth() + + +def test_tth_to_q(): + # valid tth values including edge cases when tth=0 or 180 degree, + # expected q vales are sin15, sin30, sin45, sin60, sin90 + 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 = [0, 0.258819, 0.5, 0.707107, 0.866025, 1] + assert np.allclose(actual_q, expected_q) + + +def test_tth_to_q_bad(): + # invalid tth values > 180 degree + actual = Diffraction_object(wavelength=4 * np.pi) + setattr(actual, "on_tth", [[0, 30, 60, 90, 120, 181], [1, 2, 3, 4, 5, 6]]) + with pytest.raises(ValueError): + actual.tth_to_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(): + # valid tth values including edge cases when tth=0 or 180 degree, + # expected d values are DMAX=100, 1/sin15, 1/sin30, 1/sin45, 1/sin60, 1/sin90, in reverse order + 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 = [1, 1.1547, 1.41421, 2, 3.8637, 100] + assert np.allclose(actual_d, expected_d) + + +def test_tth_to_d_bad(): + # invalid tth values > 180 degree + actual = Diffraction_object(wavelength=2) + setattr(actual, "on_tth", [[0, 30, 60, 90, 120, 181], [1, 2, 3, 4, 5, 6]]) + with pytest.raises(ValueError): + actual.tth_to_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) + + +def test_d_to_tth_bad(): + # invalid d values when arcsin value is not in the range of [-1, 1] + actual = Diffraction_object(wavelength=2) + setattr(actual, "on_d", [[0, 10, 20, 30], [1, 2, 3, 4]]) + with pytest.raises(ValueError): + actual.d_to_tth() + + +# Inputs cover a good valid range and provide the same expected values, ensuring conversion precision. +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)