custom scattering plots

Additional updates for plotting and display of custom scattering factors.

classes_crystal.py
 - Change Atoms.type to always be array type with <U8

classes_scattering.py
 - update scattering_factors

classes_plotting.py
 - added plot_scattering_factors

All tests pass.

Author: DanPorter

Dans_Diffraction/classes_crystal.py

@@ -51,6 +51,7 @@
 22/05/23 3.2.4  Added Symmetry.wyckoff_label(), Symmetry.spacegroup_dict
 06/05/24 3.3.0  Symmetry.from_cif now loads operations from find_spacegroup if not already loaded
 06/04/25 3.3.1  scale parameter of superlattice improved
+15/09/25 3.3.2  Atoms.type changed to always be array type
 
 @author: DGPorter
 """
@@ -1088,6 +1089,7 @@ class Atoms:
         "_atom_site_fract_y",
         "_atom_site_fract_z",
     ]
+    _type_str_fmt = '<U8'
 
     def __init__(self, u=[0], v=[0], w=[0], type=None,
                  label=None, occupancy=None, uiso=None, mxmymz=None):
@@ -1100,9 +1102,9 @@ def __init__(self, u=[0], v=[0], w=[0], type=None,
         # ---Defaults---
         # type
         if type is None:
-            self.type = np.asarray([self._default_atom] * Natoms)
+            self.type = np.asarray([self._default_atom] * Natoms, dtype=self._type_str_fmt)
         else:
-            self.type = np.asarray(type, dtype=str).reshape(-1)
+            self.type = np.asarray(type, dtype=self._type_str_fmt).reshape(-1)
         # label
         if label is None:
             self.label = self.type.copy()
@@ -1136,7 +1138,7 @@ def __call__(self, u=[0], v=[0], w=[0], type=None,
 
     def __getitem__(self, idx):
         if isinstance(idx, str):
-            idx = self.label.index(idx)
+            idx = list(self.label).index(idx)
         return self.atom(idx)
 
     def fromcif(self, cifvals):
@@ -1221,8 +1223,8 @@ def fromcif(self, cifvals):
         self.u = u
         self.v = v
         self.w = w
-        self.type = element
-        self.label = label
+        self.type = np.array(element, dtype=self._type_str_fmt)
+        self.label = np.array(label, dtype=self._type_str_fmt)
         self.occupancy = occ
         self.uiso = uiso
         self.mx = mx
@@ -1276,19 +1278,21 @@ def atom(self, idx):
             return atoms[0]
         return atoms
 
-    def changeatom(self, idx=None, u=None, v=None, w=None, type=None,
+    def changeatom(self, idx, u=None, v=None, w=None, type=None,
                    label=None, occupancy=None, uiso=None, mxmymz=None):
         """
-        Change an atoms properties
-        :param idx:
-        :param u:
-        :param v:
-        :param w:
-        :param type:
-        :param label:
-        :param occupancy:
-        :param uiso:
-        :param mxmymz:
+        Change an atom site's properties.
+        If properties are given as None, they are not changed.
+
+        :param idx: atom array index
+        :param u: atomic position u in relative coordinates along basis vector a
+        :param v: atomic position u in relative coordinates along basis vector b
+        :param w: atomic position u in relative coordinates along basis vector c
+        :param type: atomic element type
+        :param label: atomic site label
+        :param occupancy: atom site occupancy
+        :param uiso: atom site isotropic thermal parameter
+        :param mxmymz: atom site magnetic vector (mx, my, mz)
         :return: None
         """
 
@@ -1307,7 +1311,7 @@ def changeatom(self, idx=None, u=None, v=None, w=None, type=None,
         if label is not None:
             old_labels = list(self.label)
             old_labels[idx] = label
-            self.label = np.array(old_labels)
+            self.label = np.array(old_labels, dtype=self._type_str_fmt)
 
         if occupancy is not None:
             self.occupancy[idx] = occupancy
@@ -1433,11 +1437,11 @@ def remove_duplicates(self, min_distance=0.01, all_types=False):
         """
 
         uvw = self.uvw()
-        type = self.type
+        atom_type = self.type
         rem_atom_idx = []
-        for n in range(0, len(type) - 1):
+        for n in range(0, len(atom_type) - 1):
             match_position = fg.mag(uvw[n, :] - uvw[n + 1:, :]) < min_distance
-            match_type = type[n] == type[n + 1:]
+            match_type = atom_type[n] == atom_type[n + 1:]
             if all_types:
                 rem_atom_idx += list(1 + n + np.where(match_position)[0])
             else:

Dans_Diffraction/classes_plotting.py

@@ -1378,6 +1378,45 @@ def plot_ms_azimuth(self, hkl, energy_kev, azir=[0, 0, 1], pv=[1, 0], numsteps=3
         fp.labels(ttl, r'$\psi$ (deg)', 'Intensity')
         #plt.subplots_adjust(bottom=0.2)
 
+    def plot_scattering_factors(self, q_max=4, energy_range=None, q_range=None):
+        """
+        Plot atomic scattering factors across wavevector or energy
+
+        if q_range has more values than energy_range, figures will plot scattering factor vs Q
+        for each energy.
+        if energy_range has more values than q_range, figures will plot scattering factor vs energy
+        for each value of Q.
+
+        :param q_max: use a q_range of 0-q_max
+        :param energy_range: energy range in keV
+        :param q_range: range of wavecectors, or None to use q_max
+        :return: None
+        """
+
+        if q_range is None:
+            q_range = np.arange(0, q_max, 0.01)
+        atom_types, atom_idx = np.unique(self.xtl.Structure.type, return_index=True)
+        atom_scattering_factors = self.xtl.Scatter.scattering_factors(qmag=q_range, energy_kev=energy_range)
+        ttl = f"{self.xtl.Scatter._scattering_type}"
+        # plot multiple figures of lower dimension
+        if atom_scattering_factors.shape[2] > atom_scattering_factors.shape[0]: # energy_range > q_range
+            # different figures for different values of Q
+            for n in range(atom_scattering_factors.shape[0]):
+                plt.figure(figsize=self._figure_size, dpi=self._figure_dpi)
+                nttl = ttl + f" Q={q_range[n]:.2g}" + r" $\AA^{-1}$"
+                for atom, idx in zip(atom_types, atom_idx):
+                    plt.plot(energy_range, atom_scattering_factors[n, idx, :], label=atom)
+                fp.labels(nttl, 'Energy [keV]', 'Atomic scattering factor', legend=True)
+        else:
+            # different figures for different values of E
+            for n in range(atom_scattering_factors.shape[2]):
+                plt.figure(figsize=self._figure_size, dpi=self._figure_dpi)
+                nttl = ttl + (f" E = {energy_range[n]} keV" if atom_scattering_factors.shape[2] > 1 else "")
+                for atom, idx in zip(atom_types, atom_idx):
+                    plt.plot(q_range, np.abs(atom_scattering_factors[:, idx, n]), label=atom)
+                fp.labels(nttl, r'Q $\AA^{-1}$', 'Atomic scattering factor', legend=True)
+
+
     r''' Remove tensor_scattering 26/05/20
     def tensor_scattering_azimuth(self, atom_label, hkl, energy_kev, azir=[0, 0, 1], process='E1E1',
                                   rank=2, time=+1, parity=+1, mk=None, lk=None, sk=None):

Dans_Diffraction/classes_scattering.py

@@ -41,6 +41,7 @@
 16/05/24 2.3.4  Added printed progress bar to generate_intensity_cut during convolusion
 17/05/24 2.3.4  Changed generate_intensity_cut to make it much faster
 23/12/24 2.3.5  Added polarised neutron options
+15/09/25 2.4.0  Added custom scattering factors
 
 @author: DGPorter
 """
@@ -54,7 +55,7 @@
 from . import multiple_scattering as ms
 # from . import tensor_scattering as ts  # Removed V1.7
 
-__version__ = '2.3.5'
+__version__ = '2.4.0'
 
 
 class Scattering:
@@ -1635,13 +1636,14 @@ def scatteringbasis(self, hkl, azim_zero=(1, 0, 0), psi=0):
         Jhat_psi = fg.norm(np.cross(Qhat, Ihat_psi))
         return np.vstack([Ihat_psi, Jhat_psi, Qhat])
 
-    def scattering_factors(self, hkl, energy_kev=None):
+    def scattering_factors(self, hkl=(0, 0, 0), energy_kev=None, qmag=None):
         """Return the scattering factors[n, m] for each reflection, hkl[n,3] and each atom [m]"""
 
         if energy_kev is None:
             energy_kev = self._energy_kev
 
-        qmag = self.xtl.Cell.Qmag(hkl)
+        if qmag is None:
+            qmag = self.xtl.Cell.Qmag(hkl)
         # Scattering factors
         ff = fs.scattering_factors(
             scattering_type=self._scattering_type,
@@ -1733,23 +1735,25 @@ def print_scattering_factor_coefficients(self):
         else:
             raise Exception(f"Unknown scattering type: {scattering_type}")
 
-        coefs = fc.scattering_factor_coefficients(*self.xtl.Structure.type, table=table)
-        if not np.any(coefs[:, -1]):
-            coefs = coefs[:, :-1]  # trim final column if all zeros
-        n_doubles = coefs.shape[1] // 2
-        n_singles = coefs.shape[1] % 2
+        element_coefs = fc.scattering_factor_coefficients_custom(*self.xtl.Structure.type, default_table=table)
+
 
         structrue = self.xtl.Structure
         structure_list = zip(structrue.label, structrue.type, structrue.u, structrue.v, structrue.w)
         out = f"{scattering_type}\n"
         for n, (label, el, u, v, w) in enumerate(structure_list):
             out += f"{n:3}  {label:6} {el:5}: {u:8.3g}, {v:8.3g}, {w:8.3g}  : "
-            out += ', '.join(f"a{d} = {coefs[n, 2*d]:.3f}, A{d} = {coefs[n, 2*d+1]:.3f}" for d in range(n_doubles))
-            out += (', ' * int(n_doubles > 0)) + (f"b = {coefs[n, -1]:.3f}" * n_singles)
+
+            coefs = element_coefs[el]
+            if abs(coefs[-1]) < 1e-6:
+                coefs = coefs[:-1]  # trim final column if all zeros
+            n_doubles = len(coefs) // 2
+            n_singles = len(coefs) % 2
+            out += ', '.join(f"a{d} = {coefs[2*d]:.3f}, A{d} = {coefs[2*d+1]:.3f}" for d in range(n_doubles))
+            out += (', ' * int(n_doubles > 0)) + (f"b = {coefs[-1]:.3f}" * n_singles)
             out += '\n'
         return out
 
-
     def old_intensity(self, HKL, scattering_type=None):
         """
         Calculate the squared structure factor for the given HKL

Dans_Diffraction/functions_crystallography.py

@@ -1317,7 +1317,7 @@ def scattering_factor_coefficients(*elements, table='itc'):
         'sears': scattering_factor_coefficients_neutron_sears
     }
     table = table.lower()
-    if table.lower() not in tables:
+    if table not in tables:
         raise ValueError(f'Unknown scattering factor table: {table}')
     return tables[table](*elements)
 
@@ -1486,8 +1486,13 @@ def custom_scattering_factor_resonant(elements, q_mag, energy_kev, default_table
         tab_en, tab_f1, tab_f2 = tables[el]
         f0 = analytical_scattering_factor(q_mag, *coefs)  # (n_q, )
         f1 = np.interp(energy_kev, tab_en, tab_f1)  # (n_en, )
-        f2 = np.interp(energy_kev, tab_en, tab_f2)  # (n_ej, )
-        qff[:, n, :] = np.array([[_f0 + _f1 - 1j * _f2 for _f1, _f2 in zip(f1, f2)] for _f0 in f0])
+        f2 = np.interp(energy_kev, tab_en, tab_f2)  # (n_en, )
+        qff[:, n, :] = np.array([
+            [
+                _f0 if np.isnan(_f1) or np.isnan(_f2) else _f0 + (_f1 - 1j * _f2)
+                for _f1, _f2 in zip(f1, f2)  # loop over energy
+            ] for _f0 in f0  # loop over q values
+        ])
     return qff
 
 

Examples/example_custom_scattering_factors.py

@@ -21,8 +21,18 @@
     [-0.4833, -3.7933, -3.7862, -3.7791, -3.772, -3.765, -3.7579, -3.7508, -3.7437, -3.7367]  # f2
 ])
 dif.fc.add_custom_form_factor_coefs('Co3+', *coefs, dispersion_table=table)
-xtl.Structure.changeatom(3, type='Co3+')
+xtl.Atoms.changeatom(1, type='Co3+')
+xtl.generate_structure()
+xtl.Structure.type[3] = 'Co'
 
+
+# Display scattering factors
+print(xtl.Scatter.print_scattering_factor_coefficients())
+
+xtl.Scatter.setup_scatter(scattering_type='custom', output=False)
+xtl.Plot.plot_scattering_factors(q_max=6)
+
+# Calculate structure factors
 en = np.arange(7.7, 7.8, 0.01)
 ref = [0, 0, 3]
 inten1 = xtl.Scatter.xray_dispersion(ref, en)