Fixes for #33

functions_scattering.py
 - changed method of structure factor calculation in 'magnetic neutron' and 'magnetic xray' to give unpolarised structure factor

functions_crystallography.py
 - Changed input of magnetic_form_factor to *elements.
 - add new magnetic form factor functions:
   - magnetic_ff_symbol
   - magnetic_ff_coefs
   - print_magnetic_ff_coefs
   - last_df_orbital
   - magnetic_ff_g_factor
   - magnetic_ff_j2j0_ratio

classes_crystal.py
 - Add method scattering_factor_coefficients to Atoms.

classes_properties.py
 - Add scattering factors method to Properties

classes_scattering.py
 - add atom scattering factor methods
 - add use_waaskirf and use_sears to setup_scatter

test_data_tables.py
 - added test cases for magnetic form factors

test_structure_factors.py
 - added test cases for magnetic scattering

All tests pass.

Author: DanPorter

Dans_Diffraction/classes_crystal.py

@@ -1588,6 +1588,10 @@ def mass_fraction(self):
 
         return weights / total_weight
 
+    def scattering_factor_coefficients(self, table='itc'):
+        """Return scattering factor coefficients for the elements"""
+        return fc.scattering_factor_coefficients(*self.type, table=table)
+
     def info(self, idx=None, type=None):
         """
         Prints properties of all atoms

Dans_Diffraction/classes_properties.py

@@ -31,6 +31,7 @@
 
 from . import functions_general as fg
 from . import functions_crystallography as fc
+from . import functions_scattering as fs
 from .classes_orbitals import CrystalOrbitals
 
 __version__ = '2.0'
@@ -119,7 +120,30 @@ def magnetic_form_factor(self, hkl):
         :return: [nxm] array of scattering factors for each atom and reflection
         """
         qmag = self.xtl.Cell.Qmag(hkl)
-        return fc.magnetic_form_factor(self.xtl.Structure.type, qmag)
+        return fc.magnetic_form_factor(*self.xtl.Structure.type, qmag=qmag)
+
+    def scattering_factors(self, scattering_type, hkl, energy_kev=None,
+                           use_sears=False, use_wasskirf=False):
+        """
+        Return an array of scattering factors based on the radiation
+        :param scattering_type: str radiation, see "get_scattering_function()"
+        :param hkl: [mx1] or None, float array of wavevector magnitudes for reflections
+        :param energy_kev: [ox1] or None, float array of energies in keV
+        :param use_sears: if True, use neutron scattering lengths from ITC Vol. C, By V. F. Sears
+        :param use_wasskirf: if True, use x-ray scattering factors from Waasmaier and Kirfel
+        :return: [nxmxo] array of scattering factors
+        """
+        qmag = self.xtl.Cell.Qmag(hkl)
+        # Scattering factors
+        ff = fs.scattering_factors(
+            scattering_type=scattering_type,
+            atom_type=self.xtl.Structure.type,
+            qmag=qmag,
+            enval=energy_kev,
+            use_sears=use_sears,
+            use_wasskirf=use_wasskirf,
+        )
+        return np.squeeze(ff)
 
     def xray_edges(self):
         """

Dans_Diffraction/classes_scattering.py

@@ -183,6 +183,7 @@ def __str__(self):
         return out
 
     def setup_scatter(self, scattering_type=None, energy_kev=None, energy_mev=None, wavelength_a=None,
+                      use_sears=None, use_waaskirf=None,
                       powder_units=None, powder_pixels=None, powder_lorentz=None, powder_overlap=None,
                       int_hkl=None, specular=None, parallel=None, theta_offset=None,
                       min_theta=None, max_theta=None, min_twotheta=None, max_twotheta=None,
@@ -233,6 +234,12 @@ def setup_scatter(self, scattering_type=None, energy_kev=None, energy_mev=None,
             else:
                 self._energy_kev = fc.wave2energy(wavelength_a)
 
+        if use_sears is not None:
+            self._use_sears_scattering_lengths = use_sears
+
+        if use_waaskirf is not None:
+            self._use_waaskirf_scattering_factor = use_waaskirf
+
         if powder_units is not None:
             self._powder_units = powder_units
 
@@ -514,7 +521,7 @@ def structure_factor(self, hkl=None, scattering_type=None, int_hkl=None, **kwarg
 
                     # Get magnetic form factors
                     if self._use_magnetic_form_factor:
-                        mf = fc.magnetic_form_factor(atom_type, qmag)
+                        mf = fc.magnetic_form_factor(*atom_type, qmag=qmag)
                     else:
                         mf = None
 
@@ -1069,7 +1076,7 @@ def magnetic_neutron(self, HKL):
 
         # Get magnetic form factors
         if self._use_magnetic_form_factor:
-            ff = fc.magnetic_form_factor(atom_type, Qmag)
+            ff = fc.magnetic_form_factor(*atom_type, qmag=Qmag)
         else:
             ff = np.ones([len(HKL), Nat])
 
@@ -1135,7 +1142,7 @@ def xray_magnetic(self, HKL):
 
         # Get magnetic form factors
         if self._use_magnetic_form_factor:
-            ff = fc.magnetic_form_factor(atom_type, Qmag)
+            ff = fc.magnetic_form_factor(*atom_type, qmag=Qmag)
         else:
             ff = np.ones([len(HKL), Nat])
 
@@ -1359,7 +1366,7 @@ def xray_nonresonant_magnetic(self, HKL, energy_kev=None, azim_zero=(1, 0, 0), p
 
         # Get magnetic form factors
         if self._use_magnetic_form_factor:
-            ff = fc.magnetic_form_factor(atom_type, Qmag)
+            ff = fc.magnetic_form_factor(*atom_type, qmag=Qmag)
         else:
             ff = np.ones([len(HKL), len(uvw)])
 
@@ -1628,6 +1635,24 @@ 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):
+        """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)
+        # Scattering factors
+        ff = fs.scattering_factors(
+            scattering_type=self._scattering_type,
+            atom_type=self.xtl.Structure.type,
+            qmag=qmag,
+            enval=energy_kev,
+            use_sears=self._use_sears_scattering_lengths,
+            use_wasskirf=self._use_waaskirf_scattering_factor,
+        )
+        return ff
+
     def print_scattering_coordinates(self, hkl, azim_zero=(1, 0, 0), psi=0):
         """
         Transform magnetic vector into components within the scattering plane
@@ -1675,7 +1700,56 @@ def print_intensity(self, HKL):
             vals=(HKL[n][0],HKL[n][1],HKL[n][2],IN[n],IX[n],INM[n],IXM[n],IXRss[n],IXRsp[n],IXRps[n],IXRpp[n])
             outstr += fmt % vals
         return outstr
-    
+
+    def print_atom_scattering_factors(self, hkl, energy_kev=None):
+        """show scattering factors for each atom for each reflection"""
+        qmag = self.xtl.Cell.Qmag(hkl)
+        ff = self.scattering_factors(hkl, energy_kev)
+        hkl = np.asarray(hkl, dtype=float).reshape([-1, 3])
+
+        out = f"{self._scattering_type}\n"
+        for n, (h, k, l) in enumerate(hkl):
+            out += f"({h:.3g}, {k:.3g}, {l:.3g})  |Q| = {qmag[n]:.4g} A-1\n"
+            for m, ele in enumerate(self.xtl.Structure.type):
+                out += f"  {ele:5}: {ff[n, m]: 12.5f}\n"
+        return out
+
+    def print_scattering_factor_coefficients(self):
+        """generate string of scattering factor coefficients for each atom"""
+        scattering_type = fs.get_scattering_type(self._scattering_type)
+
+        if 'neutron' in scattering_type:
+            if self._use_sears_scattering_lengths:
+                table = 'sears'
+            else:
+                table = 'ndb'
+        elif 'electron' in scattering_type:
+            table = 'peng'
+        elif 'xray' in scattering_type:
+            if self._use_waaskirf_scattering_factor:
+                table = 'waaskirf'
+            else:
+                table = 'itc'
+        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
+
+        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)
+            out += '\n'
+        return out
+
+
     def old_intensity(self, HKL, scattering_type=None):
         """
         Calculate the squared structure factor for the given HKL
@@ -2303,48 +2377,50 @@ def print_atomic_contributions(self, HKL):
         """
         Prints the atomic contributions to the structure factor
         """
-        
-        HKL = np.asarray(np.rint(HKL),dtype=float).reshape([-1,3])
+
+        HKL = np.asarray(np.rint(HKL), dtype=float).reshape([-1, 3])
         Nref = len(HKL)
-        
+
         # Calculate the full intensity
         I = self.intensity(HKL)
-        
+
         # Calculate the structure factors of the symmetric atomic sites
         base_label = self.xtl.Atoms.label
-        uvw,type,label,occ,uiso,mxmymz = self.xtl.Structure.get()
-        
+        uvw, atom_type, label, occ, uiso, mxmymz = self.xtl.Structure.get()
+
         Qmag = self.xtl.Cell.Qmag(HKL)
-        
+
         # Get atomic form factors
-        ff = fc.xray_scattering_factor(type,Qmag)
-        
+        # ff = fc.xray_scattering_factor(atom_type, Qmag)
+        ff = self.scattering_factors(HKL)
+
         # Get Debye-Waller factor
-        dw = fc.debyewaller(uiso,Qmag)
-        
+        dw = fc.debyewaller(uiso, Qmag)
+
         # Calculate dot product
-        dot_KR = np.dot(HKL,uvw.T)
-        
+        dot_KR = np.dot(HKL, uvw.T)
+
         # Calculate structure factor
-        SF =  ff*dw*occ*np.exp(1j*2*np.pi*dot_KR)
-        
+        SF = ff * dw * occ * np.exp(1j * 2 * np.pi * dot_KR)
+
         # Sum structure factors of each base label in atoms
-        SFbase = np.zeros([len(HKL),len(base_label)],dtype=np.complex128)
+        SFbase = np.zeros([len(HKL), len(base_label) + 1], dtype=complex)
         for n in range(len(base_label)):
             label_idx = label == base_label[n]
-            SFbase[:,n] = np.sum(SF[:,label_idx],axis=1)
-        
+            SFbase[:, n] = np.sum(SF[:, label_idx], axis=1)
+        SFbase[:, -1] = SFbase.sum(axis=1)  # add the sum
+
         # Get the real part of the structure factor
-        #SFtot = np.sqrt(np.real(SF * np.conj(SF)))
+        # SFtot = np.sqrt(np.real(SF * np.conj(SF)))
         SFrel = np.real(SFbase)
         SFimg = np.imag(SFbase)
-        
+
         # Generate the results
         outstr = ''
-        outstr+= '( h, k, l) Intensity' + ' '.join(['%12s    '%x for x in base_label])+'\n'
+        outstr += '( h, k, l) Intensity ' + ' '.join(['%12s    ' % x for x in base_label]) + '      Total SF\n'
         for n in range(Nref):
-            ss = ' '.join(['%6.1f + i%-6.1f' % (x,y) for x,y in zip(SFrel[n],SFimg[n])])
-            outstr+= '(%2.0f,%2.0f,%2.0f) %9.2f    %s\n' % (HKL[n,0],HKL[n,1],HKL[n,2],I[n],ss)
+            ss = ' '.join(['%6.1f + i%-6.1f' % (x, y) for x, y in zip(SFrel[n], SFimg[n])])
+            outstr += '(%2.0f,%2.0f,%2.0f) %9.2f   %s\n' % (HKL[n, 0], HKL[n, 1], HKL[n, 2], I[n], ss)
         return outstr
 
     def print_symmetry_contributions(self, HKL):
@@ -2367,7 +2443,8 @@ def print_symmetry_contributions(self, HKL):
         # Calculate the structure factors
         uvw, type, label, occ, uiso, mxmymz = self.xtl.Structure.get()
         Qmag = self.xtl.Cell.Qmag(HKL)
-        ff = fc.xray_scattering_factor(type, Qmag)
+        # ff = fc.xray_scattering_factor(type, Qmag)
+        ff = self.scattering_factors(HKL)
         dw = fc.debyewaller(uiso, Qmag)
         dot_KR = np.dot(HKL, uvw.T)
         phase = np.exp(1j * 2 * np.pi * dot_KR)