Merge pull request #102 from isabelwood100/Izzy_diffsims_diffgen

bugfix for calculate_profile_data

diffsims/generators/diffraction_generator.py

@@ -35,6 +35,7 @@
     get_points_in_sphere,
     get_vectorized_list_for_atomic_scattering_factors,
     is_lattice_hexagonal,
+    get_intensities_params,
 )
 from diffsims.utils.fourier_transform import from_recip
 
@@ -142,34 +143,37 @@ def calculate_ed_data(
         r_sphere = 1 / wavelength
         r_spot = np.sqrt(np.sum(np.square(cartesian_coordinates[:, :2]), axis=1))
         z_sphere = -np.sqrt(r_sphere ** 2 - r_spot ** 2) + r_sphere
-        proximity = np.absolute(z_sphere - cartesian_coordinates[:, 2])
-        intersection = proximity < max_excitation_error
+        excitation_error = np.absolute(z_sphere - cartesian_coordinates[:, 2])
+        intersection = excitation_error < max_excitation_error
         # Mask parameters corresponding to excited reflections.
         intersection_coordinates = cartesian_coordinates[intersection]
-        intersection_indices = spot_indices[intersection]
-        proximity = proximity[intersection]
+        g_indices = spot_indices[intersection]
+        excitation_error = excitation_error[intersection]
         g_hkls = spot_distances[intersection]
+        multiplicites = np.ones_like(g_hkls)
+
+        shape_factor = 1 - (excitation_error / max_excitation_error)
 
         # Calculate diffracted intensities based on a kinematical model.
         intensities = get_kinematical_intensities(
             structure,
-            intersection_indices,
+            g_indices,
             g_hkls,
-            proximity,
-            max_excitation_error,
             debye_waller_factors,
+            multiplicites,
             scattering_params,
+            shape_factor,
         )
 
         # Threshold peaks included in simulation based on minimum intensity.
         peak_mask = intensities > 1e-20
         intensities = intensities[peak_mask]
         intersection_coordinates = intersection_coordinates[peak_mask]
-        intersection_indices = intersection_indices[peak_mask]
+        g_indices = g_indices[peak_mask]
 
         return DiffractionSimulation(
             coordinates=intersection_coordinates,
-            indices=intersection_indices,
+            indices=g_indices,
             intensities=intensities,
             with_direct_beam=with_direct_beam,
         )
@@ -211,79 +215,54 @@ def calculate_profile_data(
         latt = structure.lattice
         is_hex = is_lattice_hexagonal(latt)
 
-        (
-            coeffs,
-            fcoords,
-            occus,
-            dwfactors,
-        ) = get_vectorized_list_for_atomic_scattering_factors(
-            structure, {}, scattering_params=scattering_params
-        )
-
         # Obtain crystallographic reciprocal lattice points within range
         recip_latt = latt.reciprocal()
         spot_indices, _, spot_distances = get_points_in_sphere(
             recip_latt, reciprocal_radius
         )
 
-        peaks = {}
-        mask = np.logical_not((np.any(spot_indices, axis=1) == 0))
-
-        for hkl, g_hkl in zip(spot_indices[mask], spot_distances[mask]):
-            # Force miller indices to be integers.
-            hkl = [int(round(i)) for i in hkl]
-
-            d_hkl = 1 / g_hkl
-
-            # Bragg condition
-            # theta = asin(wavelength * g_hkl / 2)
-
-            # s = sin(theta) / wavelength = 1 / 2d = |ghkl| / 2 (d =
-            # 1/|ghkl|)
-            s = g_hkl / 2
-
-            # Store s^2 since we are using it a few times.
-            s2 = s ** 2
-
-            # Vectorized computation of g.r for all fractional coords and
-            # hkl.
-            g_dot_r = np.dot(fcoords, np.transpose([hkl])).T[0]
-
-            # Highly vectorized computation of atomic scattering factors.
-            fs = np.sum(coeffs[:, :, 0] * np.exp(-coeffs[:, :, 1] * s2), axis=1)
-
-            dw_correction = np.exp(-dwfactors * s2)
-
-            # Structure factor = sum of atomic scattering factors (with
-            # position factor exp(2j * pi * g.r and occupancies).
-            # Vectorized computation.
-            f_hkl = np.sum(fs * occus * np.exp(2j * np.pi * g_dot_r) * dw_correction)
+        ##spot_indicies is a numpy.array of the hkls allowd in the recip radius
+        unique_hkls, multiplicites, g_hkls = get_intensities_params(
+            recip_latt, reciprocal_radius
+        )
+        g_indices = unique_hkls
+        debye_waller_factors = self.debye_waller_factors
+        excitation_error = None
+        max_excitation_error = None
+        g_hkls_array = np.asarray(g_hkls)
 
-            # Intensity for hkl is modulus square of structure factor.
-            i_hkl = (f_hkl * f_hkl.conjugate()).real
+        i_hkl = get_kinematical_intensities(
+            structure,
+            g_indices,
+            g_hkls_array,
+            debye_waller_factors,
+            multiplicites,
+            scattering_params,
+            shape_factor=1,
+        )
 
-            # two_theta = degrees(2 * theta)
+        if is_hex:
+            # Use Miller-Bravais indices for hexagonal lattices.
+            g_indices = (g_indices[0], g_indices[1], -g_indices[0] - g_indices[1], g_indices[2])
 
-            if is_hex:
-                # Use Miller-Bravais indices for hexagonal lattices.
-                hkl = (hkl[0], hkl[1], -hkl[0] - hkl[1], hkl[2])
+        hkls_labels = ["".join([str(int(x)) for x in xs]) for xs in unique_hkls]
 
-            peaks[g_hkl] = [i_hkl, [tuple(hkl)], d_hkl]
+        peaks = {}
+        for l, i, g in zip(hkls_labels, i_hkl, g_hkls):
+            peaks[l] = [i, g]
 
         # Scale intensities so that the max intensity is 100.
+
         max_intensity = max([v[0] for v in peaks.values()])
         x = []
         y = []
         hkls = []
-        d_hkls = []
-        for k in sorted(peaks.keys()):
+        for k in peaks.keys():
             v = peaks[k]
-            fam = get_unique_families(v[1])
-            if v[0] / max_intensity * 100 > minimum_intensity:
-                x.append(k)
+            if v[0] / max_intensity * 100 > minimum_intensity and (k != "000"):
+                x.append(v[1])
                 y.append(v[0])
-                hkls.append(fam)
-                d_hkls.append(v[2])
+                hkls.append(k)
 
         y = y / max(y) * 100
 

diffsims/sims/diffraction_simulation.py

@@ -20,6 +20,8 @@
 import numpy as np
 from scipy import ndimage as ndi
 
+from diffsims.utils.sim_utils import get_unique_families
+
 
 class DiffractionSimulation:
     """Holds the result of a kinematic diffraction pattern simulation.
@@ -175,6 +177,9 @@ class ProfileSimulation:
         Magnitudes of scattering vectors.
     intensities : array-like, shape [n_peaks, 1]
         The kinematic intensity of the diffraction peaks.
+
+    Returns
+    -------
     hkls: [{(h, k, l): mult}] {(h, k, l): mult} is a dict of Miller
         indices for all diffracted lattice facets contributing to each
         intensity.
@@ -188,8 +193,12 @@ def __init__(self, magnitudes, intensities, hkls):
         self.intensities = intensities
         self.hkls = hkls
 
+
     def get_plot(self, g_max, annotate_peaks=True, with_labels=True, fontsize=12):
-        """Plots the diffraction profile simulation.
+
+        """Plots the diffraction profile simulation for the
+           calculate_profile_data method in DiffractionGenerator.
+
         Parameters
         ----------
         g_max : float
@@ -201,13 +210,13 @@ def get_plot(self, g_max, annotate_peaks=True, with_labels=True, fontsize=12):
         fontsize : integer
             Fontsize for peak labels.
         """
+
         ax = plt.gca()
         for g, i, hkls in zip(self.magnitudes, self.intensities, self.hkls):
-            if g <= g_max:
-                label = ", ".join([str(hkl) for hkl in hkls.keys()])
-                ax.plot([g, g], [0, i], color="k", linewidth=3, label=label)
-                if annotate_peaks:
-                    ax.annotate(label, xy=[g, i], xytext=[g, i], fontsize=fontsize)
+            label = hkls
+            ax.plot([g, g], [0, i], color="k", linewidth=3, label=label)
+            if annotate_peaks:
+                ax.annotate(label, xy=[g, i], xytext=[g, i], fontsize=fontsize)
 
             if with_labels:
                 ax.set_xlabel("A ($^{-1}$)")

diffsims/tests/test_sims/test_diffraction_simulation.py

@@ -66,22 +66,21 @@ def profile_simulation():
             ]
         ),
         hkls=[
-            {(1, 1, 1): 8},
-            {(2, 2, 0): 12},
-            {(3, 1, 1): 24},
-            {(4, 0, 0): 6},
-            {(3, 3, 1): 24},
-            {(4, 2, 2): 24},
-            {(3, 3, 3): 8, (5, 1, 1): 24},
-            {(4, 4, 0): 12},
-            {(5, 3, 1): 48},
-            {(6, 2, 0): 24},
-            {(5, 3, 3): 24},
-            {(4, 4, 4): 8},
+            (1, 1, 1),
+            (2, 2, 0),
+            (3, 1, 1),
+            (4, 0, 0),
+            (3, 3, 1),
+            (4, 2, 2),
+            (3, 3, 3),
+            (4, 4, 0),
+            (5, 3, 1),
+            (6, 2, 0),
+            (5, 3, 3),
+            (4, 4, 4),
         ],
     )
 
-
 def test_plot_profile_simulation(profile_simulation):
     profile_simulation.get_plot(g_max=1)
 

diffsims/tests/test_utils/test_sim_utils.py

@@ -19,7 +19,7 @@
 import pytest
 import numpy as np
 import diffpy
-import scipy.constants as sc
+from transforms3d.euler import euler2mat
 
 
 from diffsims.utils.sim_utils import (
@@ -42,6 +42,10 @@
     et_to_beta,
     acceleration_voltage_to_wavelength,
     diffraction_scattering_angle,
+    get_intensities_params,
+)
+from diffsims.tests.test_generators.test_diffraction_generator import (
+    make_structure
 )
 
 
@@ -145,7 +149,7 @@ def test_uvtw_to_uvw(uvtw, uvw):
     val = uvtw_to_uvw(uvtw)
     np.testing.assert_almost_equal(val, uvw)
 
-
+    
 class TestHolzCalibration:
     def test_get_holz_angle(self):
         wavelength = 2.51 / 1000
@@ -308,3 +312,29 @@ def test_array_like(self):
         assert len(scattering_angle) == 2
         sa_known = np.array([9.84e-3, 6.56e-3])
         assert np.allclose(sa_known, scattering_angle, rtol=0.001)
+
+def test_get_intensities_params(default_structure):
+    latt = default_structure.lattice
+    reciprocal_lattice = latt.reciprocal()
+    reciprocal_radius = 0.2
+    unique_hkls, multiplicites, g_hkls = get_intensities_params(reciprocal_lattice, reciprocal_radius)
+    np.testing.assert_equal(multiplicites, ([1.]))
+    np.testing.assert_equal(g_hkls, [0.0])
+    np.testing.assert_array_equal(unique_hkls, [[-0., -0.,  0.]])
+
+def test_get_kinematical_intensities(default_structure):
+    latt = default_structure.lattice
+    reciprocal_lattice = latt.reciprocal()
+    reciprocal_radius = 0.2
+    unique_hkls, multiplicites, g_hkls = get_intensities_params(reciprocal_lattice, reciprocal_radius)
+    g_hkls_array = np.asarray(g_hkls)
+    i_hkls = get_kinematical_intensities(
+                    default_structure,
+                    g_indices=unique_hkls,
+                    g_hkls_array=g_hkls_array,
+                    debye_waller_factors={1:1},
+                    multiplicites=multiplicites,
+                    scattering_params="lobato",
+                    shape_factor=1,
+                    )
+    np.testing.assert_array_almost_equal(i_hkls, ([43.0979]), decimal=4)