From 49057060fe3121d244778a6f5fe0f127d956652c Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Wed, 14 Apr 2021 16:56:59 +0100
Subject: [PATCH 01/14] Approximation=True precession simulation bugfix

---
 diffsims/generators/diffraction_generator.py | 50 ++++++++++----------
 1 file changed, 24 insertions(+), 26 deletions(-)

diff --git a/diffsims/generators/diffraction_generator.py b/diffsims/generators/diffraction_generator.py
index 8ec4ad62..34e20467 100644
--- a/diffsims/generators/diffraction_generator.py
+++ b/diffsims/generators/diffraction_generator.py
@@ -178,20 +178,20 @@ class DiffractionGenerator(object):
         The accelerating voltage of the microscope in kV.
     scattering_params : str
         "lobato" or "xtables"
-    minimum_intensity : float
-        Minimum intensity for a peak to be considered visible in the pattern
     precession_angle : float
         Angle about which the beam is precessed. Default is no precession.
-    approximate_precession : boolean
-        When using precession, whether to precisely calculate average
-        excitation errors and intensities or use an approximation. See notes.
     shape_factor_model : function or string
         A function that takes excitation_error and
         `max_excitation_error` (and potentially kwargs) and returns
         an intensity scaling factor. If None defaults to
         `shape_factor_models.linear`. A number of pre-programmed functions
         are available via strings.
-    kwargs
+    approximate_precession : boolean
+        When using precession, whether to precisely calculate average
+        excitation errors and intensities or use an approximation. See notes.
+    minimum_intensity : float
+        Minimum intensity for a peak to be considered visible in the pattern (fractional from the maximum)
+    kwargs :
         Keyword arguments passed to `shape_factor_model`.
 
     Notes
@@ -303,25 +303,23 @@ def calculate_ed_data(
         r_spot = np.sqrt(np.sum(np.square(cartesian_coordinates[:, :2]), axis=1))
         z_spot = cartesian_coordinates[:, 2]
 
-        if self.precession_angle > 0 and not self.approximate_precession:
-            # We find the average excitation error - this step can be
-            # quite expensive
-            excitation_error = _average_excitation_error_precession(
-                z_spot,
-                r_spot,
-                wavelength,
-                self.precession_angle,
-            )
+        z_sphere = -np.sqrt(r_sphere ** 2 - r_spot ** 2) + r_sphere
+        excitation_error = z_sphere - z_spot
+
+        if self.precession_angle == 0:
+            # Mask parameters corresponding to excited reflections.
+            intersection = np.abs(excitation_error) < max_excitation_error
+            intersection_coordinates = cartesian_coordinates[intersection]
+            excitation_error = excitation_error[intersection]
+            r_spot = r_spot[intersection]
+            g_indices = spot_indices[intersection]
+            g_hkls = spot_distances[intersection]
+
         else:
-            z_sphere = -np.sqrt(r_sphere ** 2 - r_spot ** 2) + r_sphere
-            excitation_error = z_sphere - z_spot
-        # Mask parameters corresponding to excited reflections.
-        intersection = np.abs(excitation_error) < max_excitation_error
-        intersection_coordinates = cartesian_coordinates[intersection]
-        excitation_error = excitation_error[intersection]
-        r_spot = r_spot[intersection]
-        g_indices = spot_indices[intersection]
-        g_hkls = spot_distances[intersection]
+            intersection_coordinates = cartesian_coordinates
+            g_indices = spot_indices
+            g_hkls = spot_distances
+
         # take into consideration rel-rods
         if self.precession_angle > 0 and not self.approximate_precession:
             shape_factor = _shape_factor_precession(
@@ -338,7 +336,7 @@ def calculate_ed_data(
                 excitation_error,
                 max_excitation_error,
                 r_spot,
-                self.precession_angle,
+                np.deg2rad(self.precession_angle),
             )
         else:
             shape_factor = self.shape_factor_model(
@@ -356,7 +354,7 @@ def calculate_ed_data(
         )
 
         # Threshold peaks included in simulation based on minimum intensity.
-        peak_mask = intensities > self.minimum_intensity
+        peak_mask = intensities > np.max(intensities) * self.minimum_intensity
         intensities = intensities[peak_mask]
         intersection_coordinates = intersection_coordinates[peak_mask]
         g_indices = g_indices[peak_mask]

From e9638af2752681e2f0199d8c397049d9f19ce126 Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Wed, 14 Apr 2021 16:58:05 +0100
Subject: [PATCH 02/14] internal renaming for readability

---
 diffsims/generators/diffraction_generator.py | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/diffsims/generators/diffraction_generator.py b/diffsims/generators/diffraction_generator.py
index 34e20467..f09d4921 100644
--- a/diffsims/generators/diffraction_generator.py
+++ b/diffsims/generators/diffraction_generator.py
@@ -285,7 +285,7 @@ def calculate_ed_data(
         # Obtain crystallographic reciprocal lattice points within `reciprocal_radius` and
         # g-vector magnitudes for intensity calculations.
         recip_latt = latt.reciprocal()
-        spot_indices, cartesian_coordinates, spot_distances = get_points_in_sphere(
+        g_indices, intersection_coordinates, g_distances = get_points_in_sphere(
             recip_latt, reciprocal_radius
         )
 
@@ -315,11 +315,6 @@ def calculate_ed_data(
             g_indices = spot_indices[intersection]
             g_hkls = spot_distances[intersection]
 
-        else:
-            intersection_coordinates = cartesian_coordinates
-            g_indices = spot_indices
-            g_hkls = spot_distances
-
         # take into consideration rel-rods
         if self.precession_angle > 0 and not self.approximate_precession:
             shape_factor = _shape_factor_precession(

From cbdd7253cb2ed64b9e2c5b36fb25a5e592cfd9d1 Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Wed, 14 Apr 2021 17:03:30 +0100
Subject: [PATCH 03/14] Further internal tidying of calculate_ed_data

---
 diffsims/generators/diffraction_generator.py | 50 ++++++++++----------
 1 file changed, 24 insertions(+), 26 deletions(-)

diff --git a/diffsims/generators/diffraction_generator.py b/diffsims/generators/diffraction_generator.py
index f09d4921..261a7712 100644
--- a/diffsims/generators/diffraction_generator.py
+++ b/diffsims/generators/diffraction_generator.py
@@ -297,8 +297,7 @@ def calculate_ed_data(
         R = euler2mat(ai, aj, ak, axes="rzxz")
         cartesian_coordinates = np.matmul(R, cartesian_coordinates.T).T
 
-        # Identify points intersecting the Ewald sphere within maximum
-        # excitation error and store the magnitude of their excitation error.
+        # Identify the excitation errors of candidate points
         r_sphere = 1 / wavelength
         r_spot = np.sqrt(np.sum(np.square(cartesian_coordinates[:, :2]), axis=1))
         z_spot = cartesian_coordinates[:, 2]
@@ -309,35 +308,34 @@ def calculate_ed_data(
         if self.precession_angle == 0:
             # Mask parameters corresponding to excited reflections.
             intersection = np.abs(excitation_error) < max_excitation_error
-            intersection_coordinates = cartesian_coordinates[intersection]
+            intersection_coordinates = intersection_coordinates[intersection]
             excitation_error = excitation_error[intersection]
             r_spot = r_spot[intersection]
-            g_indices = spot_indices[intersection]
-            g_hkls = spot_distances[intersection]
-
-        # take into consideration rel-rods
-        if self.precession_angle > 0 and not self.approximate_precession:
-            shape_factor = _shape_factor_precession(
-                intersection_coordinates[:, 2],
-                r_spot,
-                wavelength,
-                self.precession_angle,
-                self.shape_factor_model,
-                max_excitation_error,
-                **self.shape_factor_kwargs,
-            )
-        elif self.precession_angle > 0 and self.approximate_precession:
-            shape_factor = lorentzian_precession(
-                excitation_error,
-                max_excitation_error,
-                r_spot,
-                np.deg2rad(self.precession_angle),
-            )
-        else:
+            g_indices = g_indices[intersection]
+            g_hkls = g_distances[intersection]
+
+            # calculate shape factor
             shape_factor = self.shape_factor_model(
                 excitation_error, max_excitation_error, **self.shape_factor_kwargs
             )
-
+        else:
+            if self.approximate_precession:
+                shape_factor = lorentzian_precession(
+                    excitation_error,
+                    max_excitation_error,
+                    r_spot,
+                    np.deg2rad(self.precession_angle),
+                )
+            else:
+                shape_factor = _shape_factor_precession(
+                    intersection_coordinates[:, 2],
+                    r_spot,
+                    wavelength,
+                    self.precession_angle,
+                    self.shape_factor_model,
+                    max_excitation_error,
+                    **self.shape_factor_kwargs,
+                )
         # Calculate diffracted intensities based on a kinematical model.
         intensities = get_kinematical_intensities(
             structure,

From 0078663192c2402a77a020c6116682865af11733 Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Wed, 14 Apr 2021 17:06:52 +0100
Subject: [PATCH 04/14] More wording adjustments

---
 diffsims/generators/diffraction_generator.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/diffsims/generators/diffraction_generator.py b/diffsims/generators/diffraction_generator.py
index 261a7712..c47aef5b 100644
--- a/diffsims/generators/diffraction_generator.py
+++ b/diffsims/generators/diffraction_generator.py
@@ -188,7 +188,7 @@ class DiffractionGenerator(object):
         are available via strings.
     approximate_precession : boolean
         When using precession, whether to precisely calculate average
-        excitation errors and intensities or use an approximation. See notes.
+        excitation errors and intensities or use an approximation.
     minimum_intensity : float
         Minimum intensity for a peak to be considered visible in the pattern (fractional from the maximum)
     kwargs :
@@ -196,8 +196,6 @@ class DiffractionGenerator(object):
 
     Notes
     -----
-    * A full calculation is much slower and is not recommended for calculating
-      a diffraction library for precession diffraction patterns.
     * When using precession and approximate_precession=True, the shape factor
     model defaults to Lorentzian; shape_factor_model is ignored. Only with
     approximate_precession=False the custom shape_factor_model is used.

From 42bc5620d2e71a433b3c11356e921afd654da1e4 Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Wed, 14 Apr 2021 17:18:33 +0100
Subject: [PATCH 05/14] Tidying up the full integral case

---
 diffsims/generators/diffraction_generator.py | 82 +++-----------------
 1 file changed, 11 insertions(+), 71 deletions(-)

diff --git a/diffsims/generators/diffraction_generator.py b/diffsims/generators/diffraction_generator.py
index c47aef5b..7318bced 100644
--- a/diffsims/generators/diffraction_generator.py
+++ b/diffsims/generators/diffraction_generator.py
@@ -50,49 +50,8 @@
     "lorentzian": lorentzian,
 }
 
-
-def _z_sphere_precession(theta, r_spot, wavelength, phi):
-    """
-    Returns the z-coordinate in reciprocal space of the Ewald sphere at
-    distance r_spot from the origin
-
-    Parameters
-    ----------
-    theta : float
-        The azimuthal angle in degrees
-    r_spot : float
-        The projected length of the reciprocal lattice vector onto the plane
-        perpendicular to the optical axis in A^-1
-    wavelength : float
-        The electron wavelength in A
-    phi : float
-        The precession angle in degrees (angle between beam and optical axis)
-
-    Returns
-    -------
-    z : float
-        The height of the ewald sphere at the point r in A^-1
-
-    Notes
-    -----
-    * The azimuthal angle is the angle the beam is currently precessed to.
-    It is the angle between the projection of the beam and the projection of
-    the relevant diffraction spot both onto the x-y plane.
-    * In the derivation of this formula we assume that we will always integrate
-    over a full precession circle, because we do not explicitly take into
-    consideration x-y coordinates of reflections.
-    """
-    theta = np.deg2rad(theta)
-    r = 1 / wavelength
-    phi = np.deg2rad(phi)
-    return -np.sqrt(
-        r ** 2 * (1 - np.sin(phi) ** 2 * np.sin(theta) ** 2)
-        - (r_spot - r * np.sin(phi) * np.cos(theta)) ** 2
-    ) + r * np.cos(phi)
-
-
 def _shape_factor_precession(
-    z_spot, r_spot, wavelength, phi, shape_function, max_excitation, **kwargs
+    excitation_error, r_spot, wavelength, phi, shape_function, max_excitation, **kwargs
 ):
     """
     The rel-rod shape factors for reflections taking into account
@@ -100,12 +59,10 @@ def _shape_factor_precession(
 
     Parameters
     ----------
-    z_spot : np.ndarray (N,)
-        An array representing the z-coordinates of the reflections in A^-1
+    excitation_error : np.ndarray (N,)
+        An array of excitation errors
     r_spot : np.ndarray (N,)
         An array representing the distance of spots from the z-axis in A^-1
-    wavelength : float
-        The electron wavelength in A
     phi : float
         The precession angle in degrees
     shape_function : callable
@@ -130,34 +87,18 @@ def _shape_factor_precession(
     """
     shf = np.zeros(z_spot.shape)
     # loop over all spots
-    for i, (z_spot_i, r_spot_i) in enumerate(zip(z_spot, r_spot)):
+    for i, (excitation_error_i, r_spot_i) in enumerate(zip(excitation_error, r_spot)):
 
-        def integrand(phi):
-            z_sph = _z_sphere_precession(phi, r_spot_i, wavelength, phi)
-            return shape_function(z_spot_i - z_sph, max_excitation, **kwargs)
+        def integrand(theta):
+            # Equation 8 in L.Palatinus et al. Acta Cryst. (2019) B75, 512-522
+            S_zero = z_spot_i - z_sph
+            variable_term = r_spot_i*np.deg2rad(phi)*np.cos(theta)
+            return shape_function(S_zero + variable_term, max_excitation, **kwargs)
 
         # average factor integrated over the full revolution of the beam
-        shf[i] = (1 / (360)) * quad(integrand, 0, 360)[0]
+        shf[i] = (1 / (2*np.pi)) * quad(integrand, 0, 2*np.pi)[0]
     return shf
 
-
-def _average_excitation_error_precession(z_spot, r_spot, wavelength, precession_angle):
-    """
-    Calculate the average excitation error for spots
-    """
-    ext = np.zeros(z_spot.shape)
-    # loop over all spots
-    for i, (z_spot_i, r_spot_i) in enumerate(zip(z_spot, r_spot)):
-
-        def integrand(phi):
-            z_sph = _z_sphere_precession(phi, r_spot_i, wavelength, precession_angle)
-            return z_spot_i - z_sph
-
-        # average factor integrated over the full revolution of the beam
-        ext[i] = (1 / (360)) * quad(integrand, 0, 360)[0]
-    return ext
-
-
 class DiffractionGenerator(object):
     """Computes electron diffraction patterns for a crystal structure.
 
@@ -326,9 +267,8 @@ def calculate_ed_data(
                 )
             else:
                 shape_factor = _shape_factor_precession(
-                    intersection_coordinates[:, 2],
+                    excitation_error,
                     r_spot,
-                    wavelength,
                     self.precession_angle,
                     self.shape_factor_model,
                     max_excitation_error,

From 8ab94a3270691064074dd47f33e8acbb48ddf13b Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Wed, 14 Apr 2021 17:21:46 +0100
Subject: [PATCH 06/14] Revert some over zealous renaming

---
 diffsims/generators/diffraction_generator.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/diffsims/generators/diffraction_generator.py b/diffsims/generators/diffraction_generator.py
index 7318bced..377f74d6 100644
--- a/diffsims/generators/diffraction_generator.py
+++ b/diffsims/generators/diffraction_generator.py
@@ -224,7 +224,7 @@ def calculate_ed_data(
         # Obtain crystallographic reciprocal lattice points within `reciprocal_radius` and
         # g-vector magnitudes for intensity calculations.
         recip_latt = latt.reciprocal()
-        g_indices, intersection_coordinates, g_distances = get_points_in_sphere(
+        g_indices, cartesian_coordinates, g_distances = get_points_in_sphere(
             recip_latt, reciprocal_radius
         )
 
@@ -247,7 +247,7 @@ def calculate_ed_data(
         if self.precession_angle == 0:
             # Mask parameters corresponding to excited reflections.
             intersection = np.abs(excitation_error) < max_excitation_error
-            intersection_coordinates = intersection_coordinates[intersection]
+            intersection_coordinates = cartesian_coordinates[intersection]
             excitation_error = excitation_error[intersection]
             r_spot = r_spot[intersection]
             g_indices = g_indices[intersection]
@@ -258,6 +258,8 @@ def calculate_ed_data(
                 excitation_error, max_excitation_error, **self.shape_factor_kwargs
             )
         else:
+            intersection_coordinates = cartesian_coordinates #for naming simplicity
+            
             if self.approximate_precession:
                 shape_factor = lorentzian_precession(
                     excitation_error,

From bde62efee578f9d4daab098b14be890cf6548601 Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Wed, 14 Apr 2021 17:35:57 +0100
Subject: [PATCH 07/14] More tweaking for no approximation precession

---
 diffsims/generators/diffraction_generator.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/diffsims/generators/diffraction_generator.py b/diffsims/generators/diffraction_generator.py
index 377f74d6..5b55dbd0 100644
--- a/diffsims/generators/diffraction_generator.py
+++ b/diffsims/generators/diffraction_generator.py
@@ -51,7 +51,7 @@
 }
 
 def _shape_factor_precession(
-    excitation_error, r_spot, wavelength, phi, shape_function, max_excitation, **kwargs
+    excitation_error, r_spot, phi, shape_function, max_excitation, **kwargs
 ):
     """
     The rel-rod shape factors for reflections taking into account
@@ -85,13 +85,13 @@ def _shape_factor_precession(
     to the optical axis, so that the shape factor function only depends on the
     distance from each spot to the Ewald sphere parallel to the optical axis.
     """
-    shf = np.zeros(z_spot.shape)
+    shf = np.zeros(excitation_error.shape)
     # loop over all spots
     for i, (excitation_error_i, r_spot_i) in enumerate(zip(excitation_error, r_spot)):
 
         def integrand(theta):
             # Equation 8 in L.Palatinus et al. Acta Cryst. (2019) B75, 512-522
-            S_zero = z_spot_i - z_sph
+            S_zero = excitation_error_i
             variable_term = r_spot_i*np.deg2rad(phi)*np.cos(theta)
             return shape_function(S_zero + variable_term, max_excitation, **kwargs)
 
@@ -224,7 +224,7 @@ def calculate_ed_data(
         # Obtain crystallographic reciprocal lattice points within `reciprocal_radius` and
         # g-vector magnitudes for intensity calculations.
         recip_latt = latt.reciprocal()
-        g_indices, cartesian_coordinates, g_distances = get_points_in_sphere(
+        g_indices, cartesian_coordinates, g_hkls = get_points_in_sphere(
             recip_latt, reciprocal_radius
         )
 
@@ -251,7 +251,7 @@ def calculate_ed_data(
             excitation_error = excitation_error[intersection]
             r_spot = r_spot[intersection]
             g_indices = g_indices[intersection]
-            g_hkls = g_distances[intersection]
+            g_hkls = g_hkls[intersection]
 
             # calculate shape factor
             shape_factor = self.shape_factor_model(
@@ -259,7 +259,7 @@ def calculate_ed_data(
             )
         else:
             intersection_coordinates = cartesian_coordinates #for naming simplicity
-            
+
             if self.approximate_precession:
                 shape_factor = lorentzian_precession(
                     excitation_error,
@@ -271,7 +271,7 @@ def calculate_ed_data(
                 shape_factor = _shape_factor_precession(
                     excitation_error,
                     r_spot,
-                    self.precession_angle,
+                    np.deg2rad(self.precession_angle),
                     self.shape_factor_model,
                     max_excitation_error,
                     **self.shape_factor_kwargs,

From c4bb00629adfe08f399041acf3ac55bd2d1a0cf3 Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Wed, 14 Apr 2021 17:36:35 +0100
Subject: [PATCH 08/14] Fix an essental typo in precession quick sim

---
 diffsims/utils/shape_factor_models.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/diffsims/utils/shape_factor_models.py b/diffsims/utils/shape_factor_models.py
index a0bcfc7c..f9dad628 100644
--- a/diffsims/utils/shape_factor_models.py
+++ b/diffsims/utils/shape_factor_models.py
@@ -199,6 +199,6 @@ def lorentzian_precession(
     """
     sigma = np.pi / max_excitation_error
     u = sigma ** 2 * (r_spot ** 2 * precession_angle ** 2 - excitation_error ** 2) + 1
-    z = np.sqrt(u ** 2 + 4 * sigma ** 2 + excitation_error ** 2)
+    z = np.sqrt(u ** 2 + 4 * sigma ** 2 * excitation_error ** 2)
     fac = (sigma / np.pi) * np.sqrt(2 * (u + z) / z ** 2)
     return fac

From e076316caf76e1f06b0408f699c94ba449ab5931 Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Wed, 14 Apr 2021 20:59:50 +0100
Subject: [PATCH 09/14] Docstring updates

---
 diffsims/generators/diffraction_generator.py | 4 ++--
 diffsims/utils/shape_factor_models.py        | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/diffsims/generators/diffraction_generator.py b/diffsims/generators/diffraction_generator.py
index 5b55dbd0..04970f7a 100644
--- a/diffsims/generators/diffraction_generator.py
+++ b/diffsims/generators/diffraction_generator.py
@@ -64,7 +64,7 @@ def _shape_factor_precession(
     r_spot : np.ndarray (N,)
         An array representing the distance of spots from the z-axis in A^-1
     phi : float
-        The precession angle in degrees
+        The precession angle in radians
     shape_function : callable
         A function that describes the influence from the rel-rods. Should be
         in the form func(excitation_error: np.ndarray, max_excitation: float,
@@ -92,7 +92,7 @@ def _shape_factor_precession(
         def integrand(theta):
             # Equation 8 in L.Palatinus et al. Acta Cryst. (2019) B75, 512-522
             S_zero = excitation_error_i
-            variable_term = r_spot_i*np.deg2rad(phi)*np.cos(theta)
+            variable_term = r_spot_i*(phi)*np.cos(theta)
             return shape_function(S_zero + variable_term, max_excitation, **kwargs)
 
         # average factor integrated over the full revolution of the beam
diff --git a/diffsims/utils/shape_factor_models.py b/diffsims/utils/shape_factor_models.py
index f9dad628..e2a3cc19 100644
--- a/diffsims/utils/shape_factor_models.py
+++ b/diffsims/utils/shape_factor_models.py
@@ -185,7 +185,7 @@ def lorentzian_precession(
         The distance (reciprocal) from each reflection to the origin
 
     precession_angle : float
-        The beam precession angle in degrees; the angle the beam makes
+        The beam precession angle in radians; the angle the beam makes
         with the optical axis.
 
     Returns

From f5bc126641e7ad87f4abd62e4f601a8e9aa01484 Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Wed, 14 Apr 2021 21:27:18 +0100
Subject: [PATCH 10/14] update tests

---
 .../generators/test_diffraction_generator.py  | 30 ++-----------------
 1 file changed, 2 insertions(+), 28 deletions(-)

diff --git a/diffsims/tests/generators/test_diffraction_generator.py b/diffsims/tests/generators/test_diffraction_generator.py
index 5000324b..39e6ccb7 100644
--- a/diffsims/tests/generators/test_diffraction_generator.py
+++ b/diffsims/tests/generators/test_diffraction_generator.py
@@ -23,9 +23,7 @@
 from diffsims.generators.diffraction_generator import (
     DiffractionGenerator,
     AtomicDiffractionGenerator,
-    _z_sphere_precession,
     _shape_factor_precession,
-    _average_excitation_error_precession,
 )
 import diffpy.structure
 from diffsims.utils.shape_factor_models import linear, binary, sin2c, atanc, lorentzian
@@ -103,35 +101,11 @@ def probe(x, out=None, scale=None):
         v = v * abs(x[1].reshape(1, -1, 1)) < 6
         return v + 0 * x[2].reshape(1, 1, -1)
 
-
-@pytest.mark.parametrize(
-    "parameters, expected",
-    [
-        ([0, 1, 0.001, 0.5], -0.00822681491001731),
-        (
-            [0, np.array([1, 2, 20]), 0.001, 0.5],
-            np.array([-0.00822681, -0.01545354, 0.02547058]),
-        ),
-        ([180, 1, 0.001, 0.5], 0.00922693),
-    ],
-)
-def test_z_sphere_precession(parameters, expected):
-    result = _z_sphere_precession(*parameters)
-    assert np.allclose(result, expected)
-
-
 @pytest.mark.parametrize("model", [binary, linear, atanc, sin2c, lorentzian])
 def test_shape_factor_precession(model):
-    z = np.array([-0.1, 0.1])
+    excitation = np.array([-0.1, 0.1])
     r = np.array([1, 5])
-    _ = _shape_factor_precession(z, r, 0.001, 0.5, model, 0.1)
-
-
-def test_average_excitation_error_precession():
-    z = np.array([-0.1, 0.1])
-    r = np.array([1, 5])
-    _ = _average_excitation_error_precession(z, r, 0.001, 0.5)
-
+    _ = _shape_factor_precession(excitation, r, 0.5, model, 0.1)
 
 @pytest.mark.parametrize(
     "model, expected",

From 77fabc10b111dc3a588ba5e4be46d1d2f6fa4506 Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Wed, 14 Apr 2021 21:30:21 +0100
Subject: [PATCH 11/14] Changelog, and release to 0.4.2-dev

---
 CHANGELOG.md             | 4 +++-
 diffsims/release_info.py | 2 +-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9b6b3c68..4457efe5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,8 +9,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added 
 - Simulations now have a .get_as_mask() method (#154, #158)
 
-## 2021-03-15 - version 0.4.1
+### Fixed 
+- Precession simulations (#161)
 
+## 2021-03-15 - version 0.4.1
 ### Changed
 - `get_grid_beam_directions` default meshing changed to "spherified_cube_edge" from "spherified_cube_corner"
 
diff --git a/diffsims/release_info.py b/diffsims/release_info.py
index 9046656a..bd1574b3 100644
--- a/diffsims/release_info.py
+++ b/diffsims/release_info.py
@@ -1,5 +1,5 @@
 name = "diffsims"
-version = "0.5.0-dev"
+version = "0.4.2-dev"
 author = "Duncan Johnstone, Phillip Crout"
 copyright = "Copyright 2017-2021, The pyXem Developers"
 credits = [

From 91e4958f9cc37cc5c8719308fa0aa0485cbce75e Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Fri, 16 Apr 2021 15:24:20 +0100
Subject: [PATCH 12/14] Testing includes 3.9

---
 .github/workflows/build.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index c202c6d7..20e33339 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -13,7 +13,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, windows-latest, macos-latest]
-        python-version: [3.7, 3.8]
+        python-version: [3.7, 3.8,3.9]
     steps:
       - uses: actions/checkout@v2
       - name: Set up Python ${{ matrix.python-version }}
@@ -43,4 +43,4 @@ jobs:
     - name: Coveralls finished
       uses: AndreMiras/coveralls-python-action@develop
       with:
-        parallel-finished: true
\ No newline at end of file
+        parallel-finished: true

From 9382bfaf3ca2da09eca2706482a87181a9ec6252 Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Fri, 16 Apr 2021 15:25:37 +0100
Subject: [PATCH 13/14] Update CHANGELOG in prep for release

---
 CHANGELOG.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4457efe5..d1552a0e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,10 +8,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added 
 - Simulations now have a .get_as_mask() method (#154, #158)
+- Python 3.9 testing (#161)
 
 ### Fixed 
 - Precession simulations (#161)
 
+### Changed
+- Simulations now use a fractional (rather than absolute) min_intensity (#161)
+
 ## 2021-03-15 - version 0.4.1
 ### Changed
 - `get_grid_beam_directions` default meshing changed to "spherified_cube_edge" from "spherified_cube_corner"

From 054bf799411ed5da10b3cf7ed493c15b842fe4bb Mon Sep 17 00:00:00 2001
From: phillip <pc494@cam.ac.uk>
Date: Fri, 16 Apr 2021 15:29:51 +0100
Subject: [PATCH 14/14] Release version

---
 CHANGELOG.md             | 2 +-
 diffsims/release_info.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index d1552a0e..88f5adff 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,7 +4,7 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## Unreleased
+## 2021-04-16 - version 0.4.2
 
 ### Added 
 - Simulations now have a .get_as_mask() method (#154, #158)
diff --git a/diffsims/release_info.py b/diffsims/release_info.py
index bd1574b3..19c5dedf 100644
--- a/diffsims/release_info.py
+++ b/diffsims/release_info.py
@@ -1,5 +1,5 @@
 name = "diffsims"
-version = "0.4.2-dev"
+version = "0.4.2"
 author = "Duncan Johnstone, Phillip Crout"
 copyright = "Copyright 2017-2021, The pyXem Developers"
 credits = [