Source code for diffpy.utils.scattering_objects.diffraction_objects

import datetime
from copy import deepcopy

import numpy as np

from diffpy.utils.tools import get_package_info

QQUANTITIES = ["q"]
ANGLEQUANTITIES = ["angle", "tth", "twotheta", "2theta"]
DQUANTITIES = ["d", "dspace"]
XQUANTITIES = ANGLEQUANTITIES + DQUANTITIES + QQUANTITIES
XUNITS = ["degrees", "radians", "rad", "deg", "inv_angs", "inv_nm", "nm-1", "A-1"]

x_grid_emsg = (
    "objects are not on the same x-grid. You may add them using the self.add method "
    "and specifying how to handle the mismatch."
)


[docs] class Diffraction_object: def __init__(self, name="", wavelength=None): self.name = name self.wavelength = wavelength self.scat_quantity = "" self.on_q = [np.empty(0), np.empty(0)] self.on_tth = [np.empty(0), np.empty(0)] self.on_d = [np.empty(0), np.empty(0)] self._all_arrays = [self.on_q, self.on_tth] self.metadata = {} def __eq__(self, other): if not isinstance(other, Diffraction_object): return NotImplemented self_attributes = [key for key in self.__dict__ if not key.startswith("_")] other_attributes = [key for key in other.__dict__ if not key.startswith("_")] if not sorted(self_attributes) == sorted(other_attributes): return False for key in self_attributes: value = getattr(self, key) other_value = getattr(other, key) if isinstance(value, float): if ( not (value is None and other_value is None) and (value is None) or (other_value is None) or not np.isclose(value, other_value, rtol=1e-5) ): return False elif isinstance(value, list) and all(isinstance(i, np.ndarray) for i in value): if not all(np.allclose(i, j, rtol=1e-5) for i, j in zip(value, other_value)): return False else: if value != other_value: return False return True def __add__(self, other): summed = deepcopy(self) if isinstance(other, int) or isinstance(other, float) or isinstance(other, np.ndarray): summed.on_tth[1] = self.on_tth[1] + other summed.on_q[1] = self.on_q[1] + other elif not isinstance(other, Diffraction_object): raise TypeError("I only know how to sum two Diffraction_object objects") elif self.on_tth[0].all() != other.on_tth[0].all(): raise RuntimeError(x_grid_emsg) else: summed.on_tth[1] = self.on_tth[1] + other.on_tth[1] summed.on_q[1] = self.on_q[1] + other.on_q[1] return summed def __radd__(self, other): summed = deepcopy(self) if isinstance(other, int) or isinstance(other, float) or isinstance(other, np.ndarray): summed.on_tth[1] = self.on_tth[1] + other summed.on_q[1] = self.on_q[1] + other elif not isinstance(other, Diffraction_object): raise TypeError("I only know how to sum two Scattering_object objects") elif self.on_tth[0].all() != other.on_tth[0].all(): raise RuntimeError(x_grid_emsg) else: summed.on_tth[1] = self.on_tth[1] + other.on_tth[1] summed.on_q[1] = self.on_q[1] + other.on_q[1] return summed def __sub__(self, other): subtracted = deepcopy(self) if isinstance(other, int) or isinstance(other, float) or isinstance(other, np.ndarray): subtracted.on_tth[1] = self.on_tth[1] - other subtracted.on_q[1] = self.on_q[1] - other elif not isinstance(other, Diffraction_object): raise TypeError("I only know how to subtract two Scattering_object objects") elif self.on_tth[0].all() != other.on_tth[0].all(): raise RuntimeError(x_grid_emsg) else: subtracted.on_tth[1] = self.on_tth[1] - other.on_tth[1] subtracted.on_q[1] = self.on_q[1] - other.on_q[1] return subtracted def __rsub__(self, other): subtracted = deepcopy(self) if isinstance(other, int) or isinstance(other, float) or isinstance(other, np.ndarray): subtracted.on_tth[1] = other - self.on_tth[1] subtracted.on_q[1] = other - self.on_q[1] elif not isinstance(other, Diffraction_object): raise TypeError("I only know how to subtract two Scattering_object objects") elif self.on_tth[0].all() != other.on_tth[0].all(): raise RuntimeError(x_grid_emsg) else: subtracted.on_tth[1] = other.on_tth[1] - self.on_tth[1] subtracted.on_q[1] = other.on_q[1] - self.on_q[1] return subtracted def __mul__(self, other): multiplied = deepcopy(self) if isinstance(other, int) or isinstance(other, float) or isinstance(other, np.ndarray): multiplied.on_tth[1] = other * self.on_tth[1] multiplied.on_q[1] = other * self.on_q[1] elif not isinstance(other, Diffraction_object): raise TypeError("I only know how to multiply two Scattering_object objects") elif self.on_tth[0].all() != other.on_tth[0].all(): raise RuntimeError(x_grid_emsg) else: multiplied.on_tth[1] = self.on_tth[1] * other.on_tth[1] multiplied.on_q[1] = self.on_q[1] * other.on_q[1] return multiplied def __rmul__(self, other): multiplied = deepcopy(self) if isinstance(other, int) or isinstance(other, float) or isinstance(other, np.ndarray): multiplied.on_tth[1] = other * self.on_tth[1] multiplied.on_q[1] = other * self.on_q[1] elif self.on_tth[0].all() != other.on_tth[0].all(): raise RuntimeError(x_grid_emsg) else: multiplied.on_tth[1] = self.on_tth[1] * other.on_tth[1] multiplied.on_q[1] = self.on_q[1] * other.on_q[1] return multiplied def __truediv__(self, other): divided = deepcopy(self) if isinstance(other, int) or isinstance(other, float) or isinstance(other, np.ndarray): divided.on_tth[1] = other / self.on_tth[1] divided.on_q[1] = other / self.on_q[1] elif not isinstance(other, Diffraction_object): raise TypeError("I only know how to multiply two Scattering_object objects") elif self.on_tth[0].all() != other.on_tth[0].all(): raise RuntimeError(x_grid_emsg) else: divided.on_tth[1] = self.on_tth[1] / other.on_tth[1] divided.on_q[1] = self.on_q[1] / other.on_q[1] return divided def __rtruediv__(self, other): divided = deepcopy(self) if isinstance(other, int) or isinstance(other, float) or isinstance(other, np.ndarray): divided.on_tth[1] = other / self.on_tth[1] divided.on_q[1] = other / self.on_q[1] elif self.on_tth[0].all() != other.on_tth[0].all(): raise RuntimeError(x_grid_emsg) else: divided.on_tth[1] = other.on_tth[1] / self.on_tth[1] divided.on_q[1] = other.on_q[1] / self.on_q[1] return divided
[docs] def set_angles_from_list(self, angles_list): self.angles = angles_list self.n_steps = len(angles_list) - 1.0 self.begin_angle = self.angles[0] self.end_angle = self.angles[-1]
[docs] def set_qs_from_range(self, begin_q, end_q, step_size=None, n_steps=None): """ create an array of linear spaced Q-values Parameters ---------- begin_q float the beginning angle end_q float the ending angle step_size float the size of the step between points. Only specify step_size or n_steps, not both n_steps integer the number of steps. Odd numbers are preferred. Only specify step_size or n_steps, not both Returns ------- Sets self.qs self.qs array of floats the q values in the independent array """ self.qs = self._set_array_from_range(begin_q, end_q, step_size=step_size, n_steps=n_steps)
[docs] def set_angles_from_range(self, begin_angle, end_angle, step_size=None, n_steps=None): """ create an array of linear spaced angle-values Parameters ---------- begin_angle float the beginning angle end_angle float the ending angle step_size float the size of the step between points. Only specify step_size or n_steps, not both n_steps integer the number of steps. Odd numbers are preferred. Only specify step_size or n_steps, not both Returns ------- Sets self.angles self.angles array of floats the q values in the independent array """ self.angles = self._set_array_from_range(begin_angle, end_angle, step_size=step_size, n_steps=n_steps)
def _set_array_from_range(self, begin, end, step_size=None, n_steps=None): if step_size is not None and n_steps is not None: print( "WARNING: both step_size and n_steps have been given. n_steps will be used and step_size will be " "reset." ) array = np.linspace(begin, end, n_steps) elif step_size is not None: array = np.arange(begin, end, step_size) elif n_steps is not None: array = np.linspace(begin, end, n_steps) return array
[docs] def get_angle_index(self, angle): count = 0 for i, target in enumerate(self.angles): if angle == target: return i else: count += 1 if count >= len(self.angles): raise IndexError(f"WARNING: no angle {angle} found in angles list")
[docs] def insert_scattering_quantity( self, xarray, yarray, xtype, metadata={}, scat_quantity=None, name=None, wavelength=None, ): f""" insert a new scattering quantity into the scattering object Parameters ---------- xarray array-like of floats the independent variable array yarray array-like of floats the dependent variable array xtype string the type of quantity for the independent variable from {*XQUANTITIES, } metadata: dict the metadata in the form of a dictionary of user-supplied key:value pairs Returns ------- """ self.input_xtype = xtype # empty attributes have been defined in the __init__ method so only # set the attributes that are not empty to avoid emptying them by mistake if metadata: self.metadata = metadata if scat_quantity is not None: self.scat_quantity = scat_quantity if name is not None: self.name = name if wavelength is not None: self.wavelength = wavelength if xtype.lower() in QQUANTITIES: self.on_q = [np.array(xarray), np.array(yarray)] 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.set_all_arrays()
[docs] def q_to_tth(self): r""" Helper function to convert q to two-theta. By definition the relationship is: .. math:: \sin\left(\frac{2\theta}{2}\right) = \frac{\lambda q}{4 \pi} thus .. math:: 2\theta_n = 2 \arcsin\left(\frac{\lambda q}{4 \pi}\right) Parameters ---------- q : array An array of :math:`q` values wavelength : float Wavelength of the incoming x-rays Function adapted from scikit-beam. Thanks to those developers Returns ------- two_theta : array An 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))
[docs] def tth_to_q(self): r""" Helper function to convert two-theta to q By definition the relationship is .. math:: \sin\left(\frac{2\theta}{2}\right) = \frac{\lambda q}{4 \pi} thus .. math:: q = \frac{4 \pi \sin\left(\frac{2\theta}{2}\right)}{\lambda} Parameters ---------- two_theta : array An array of :math:`2\theta` values in units of degrees wavelength : float Wavelength of the incoming x-rays Function adapted from scikit-beam. Thanks to those developers. Returns ------- q : array An 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)
[docs] 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] if xtype == "tth": self.on_q[0] = self.tth_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]
def _get_original_array(self): if self.input_xtype in QQUANTITIES: return self.on_q, "q" elif self.input_xtype in ANGLEQUANTITIES: return self.on_tth, "tth" elif self.input_xtype in DQUANTITIES: return self.on_d, "d"
[docs] def scale_to(self, target_diff_object, xtype=None, xvalue=None): f""" returns a new diffraction object which is the current object but recaled in y to the target Parameters ---------- target_diff_object: Diffraction_object the diffractoin object you want to scale the current one on to xtype: string, optional. Default is Q the xtype, from {XQUANTITIES}, that you will specify a point from to scale to xvalue: float. Default is the midpoint of the array the y-value in the target at this x-value will be used as the factor to scale to. The entire array is scaled be the factor that places on on top of the other at that point. xvalue does not have to be in the x-array, the point closest to this point will be used for the scaling. Returns ------- the rescaled Diffraction_object as a new object """ scaled = deepcopy(self) if xtype is None: xtype = "q" data = self.on_xtype(xtype) target = target_diff_object.on_xtype(xtype) if xvalue is None: xvalue = data[0][0] + (data[0][-1] - data[0][0]) / 2.0 xindex = (np.abs(data[0] - xvalue)).argmin() ytarget = target[1][xindex] yself = data[1][xindex] scaled.on_tth[1] = data[1] * ytarget / yself scaled.on_q[1] = data[1] * ytarget / yself return scaled
[docs] def on_xtype(self, xtype): """ return a 2D np array with x in the first column and y in the second for x of type type Parameters ---------- xtype Returns ------- """ if xtype.lower() in ANGLEQUANTITIES: return self.on_tth elif xtype.lower() in QQUANTITIES: return self.on_q elif xtype.lower() in DQUANTITIES: return self.on_d pass
[docs] def dump(self, filepath, xtype=None): if xtype is None: xtype = " q" if xtype == "q": data_to_save = np.column_stack((self.on_q[0], self.on_q[1])) elif xtype == "tth": data_to_save = np.column_stack((self.on_tth[0], self.on_tth[1])) else: print(f"WARNING: cannot handle the xtype '{xtype}'") self.metadata.update(get_package_info("diffpy.utils", metadata=self.metadata)) self.metadata["creation_time"] = datetime.datetime.now() with open(filepath, "w") as f: f.write( f"[Diffraction_object]\nname = {self.name}\nwavelength = {self.wavelength}\n" f"scat_quantity = {self.scat_quantity}\n" ) for key, value in self.metadata.items(): f.write(f"{key} = {value}\n") f.write("\n#### start data\n") np.savetxt(f, data_to_save, delimiter=" ")