diff --git a/CHANGELOG.md b/CHANGELOG.md
index c9994e74..5666f417 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,11 +6,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 ### Changed
-- get_grid_beam_directions, now works based off of meshes
+- `get_grid_beam_directions`, now works based off of meshes
+- the arguments in the `DiffractionGenerator` constructor and the `DiffractionLibraryGenerator.get_diffraction_library` function have been shuffled so that the former captures arguments related to "the instrument/physics" while the latter captures arguments relevant to "the sample/material".
 
 ### Added
 - API reference documentation via Read The Docs: https://diffsims.rtfd.io
-- New module: "sphere_mesh_generators"
+- New module: `sphere_mesh_generators`
+- beam precession is now supported in simulating electron diffraction patterns
+- more shape factor functions have been added
 - This project now keeps a Changelog
 
 ### Removed
diff --git a/diffsims/generators/diffraction_generator.py b/diffsims/generators/diffraction_generator.py
index f9c13b10..3d09cc26 100644
--- a/diffsims/generators/diffraction_generator.py
+++ b/diffsims/generators/diffraction_generator.py
@@ -51,40 +51,48 @@
 }
 
 
-def _z_sphere_precession(phi, r_spot, wavelength, theta):
+def _z_sphere_precession(theta, r_spot, wavelength, phi):
     """
     Returns the z-coordinate in reciprocal space of the Ewald sphere at
-    distance r from the origin along the x-axis
+    distance r_spot from the origin
 
     Parameters
     ----------
-    phi : float
-        The angle the beam is currently precessed to in degrees.
-        If the beam were projected on the x-y plane, the angle
-        is the angle between this projection with the x-axis.
+    theta : float
+        The azimuthal angle in degrees
     r_spot : float
-        The distance to the point on the x-axis where we calculate z in A^-1
+        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^-1
-    theta : 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.
     """
-    phi = np.deg2rad(phi)
-    r = 1 / wavelength
     theta = np.deg2rad(theta)
+    r = 1 / wavelength
+    phi = np.deg2rad(phi)
     return -np.sqrt(
-        r ** 2 * (1 - np.sin(theta) ** 2 * np.sin(phi) ** 2)
-        - (r_spot - r * np.sin(theta) * np.cos(phi)) ** 2
-    ) + r * np.cos(theta)
+        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, precession_angle, function, max_excitation, **kwargs
+    z_spot, r_spot, wavelength, phi, shape_function, max_excitation, **kwargs
 ):
     """
     The rel-rod shape factors for reflections taking into account
@@ -98,34 +106,35 @@ def _shape_factor_precession(
         An array representing the distance of spots from the z-axis in A^-1
     wavelength : float
         The electron wavelength in A
-    precession_angle : float
+    phi : float
         The precession angle in degrees
-    function : callable
+    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,
         **kwargs)
     max_excitation : float
-        Maximum excitation angle to consider if precession_angle = 0. With
-        precession, it is rather a parameter to describe the "width" of the
-        rel-rods.
+        Parameter to describe the "extent" of the rel-rods.
 
     Other parameters
     ----------------
-    ** kwargs: passed directly to function
+    ** kwargs: passed directly to shape_function
 
     Notes
     -----
-    We calculate excitation_error as z_spot - z_sphere so that it is
+    * We calculate excitation_error as z_spot - z_sphere so that it is
     negative when the spot is outside the ewald sphere and positive when inside
     conform W&C chapter 12, section 12.6
+    * We assume that the sample is a thin infinitely wide slab perpendicular
+    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)
     # 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 function(z_spot_i - z_sph, max_excitation, **kwargs)
+            z_sph = _z_sphere_precession(phi, r_spot_i, wavelength, phi)
+            return shape_function(z_spot_i - z_sph, max_excitation, **kwargs)
 
         # average factor integrated over the full revolution of the beam
         shf[i] = (1 / (360)) * quad(integrand, 0, 360)[0]
@@ -175,7 +184,7 @@ class DiffractionGenerator(object):
         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.
+        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
@@ -190,7 +199,7 @@ class DiffractionGenerator(object):
     * 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
+    model defaults to Lorentzian; shape_factor_model is ignored. Only with
     approximate_precession=False the custom shape_factor_model is used.
     """
 
diff --git a/diffsims/utils/shape_factor_models.py b/diffsims/utils/shape_factor_models.py
index c14f4a00..fb867aa7 100644
--- a/diffsims/utils/shape_factor_models.py
+++ b/diffsims/utils/shape_factor_models.py
@@ -81,9 +81,9 @@ def sinc(excitation_error, max_excitation_error, minima_number=5):
     num = np.sin(fac * excitation_error)
     denom = fac * excitation_error
     return np.nan_to_num(
-            np.abs(np.divide(num, denom, out=np.zeros_like(num), where=denom != 0)),
-            nan=1,
-            )
+        np.abs(np.divide(num, denom, out=np.zeros_like(num), where=denom != 0)),
+        nan=1,
+    )
 
 
 def sin2c(excitation_error, max_excitation_error, minima_number=5):
@@ -105,7 +105,7 @@ def sin2c(excitation_error, max_excitation_error, minima_number=5):
     -------
     intensity : array-like or float
     """
-    return sinc(excitation_error, max_excitation_error, minima_number)**2
+    return sinc(excitation_error, max_excitation_error, minima_number) ** 2
 
 
 def atanc(excitation_error, max_excitation_error, minima_number=5):
@@ -131,15 +131,15 @@ def atanc(excitation_error, max_excitation_error, minima_number=5):
     """
     fac = np.pi * minima_number / np.abs(max_excitation_error)
     return np.nan_to_num(
-            np.arctan(fac*excitation_error)/(fac*excitation_error),
-            nan=1,
-            )
+        np.arctan(fac * excitation_error) / (fac * excitation_error),
+        nan=1,
+    )
 
 
 def lorentzian(excitation_error, max_excitation_error):
     """
     Lorentzian intensity profile that should approximate
-    the two-beam rocking curve. After [1].
+    the two-beam rocking curve. This is equation (6) in reference [1].
 
     Parameters
     ----------
@@ -160,16 +160,18 @@ def lorentzian(excitation_error, max_excitation_error):
     """
     # in the paper, sigma = pi*thickness.
     # We assume thickness = 1/max_exitation_error
-    sigma = np.pi/max_excitation_error
-    fac = sigma/(np.pi*(sigma**2*excitation_error**2+1))
+    sigma = np.pi / max_excitation_error
+    fac = sigma / (np.pi * (sigma ** 2 * excitation_error ** 2 + 1))
     return fac
 
 
-def lorentzian_precession(excitation_error, max_excitation_error,
-                          r_spot, precession_angle):
+def lorentzian_precession(
+    excitation_error, max_excitation_error, r_spot, precession_angle
+):
     """
     Intensity profile factor for a precessed beam assuming a Lorentzian
-    intensity profile for the un-precessed beam. After [1].
+    intensity profile for the un-precessed beam. This is equation (10) in
+    reference [1].
 
     Parameters
     ----------
@@ -195,8 +197,8 @@ def lorentzian_precession(excitation_error, max_excitation_error,
     ----------
     [1] L. Palatinus, P. Brázda, M. Jelínek, J. Hrdá, G. Steciuk, M. Klementová, Specifics of the data processing of precession electron diffraction tomography data and their implementation in the program PETS2.0, Acta Crystallogr. Sect. B Struct. Sci. Cryst. Eng. Mater. 75 (2019) 512–522. doi:10.1107/S2052520619007534.
     """
-    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)
-    fac = (sigma/np.pi)*np.sqrt(2*(u+z)/z**2)
+    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)
+    fac = (sigma / np.pi) * np.sqrt(2 * (u + z) / z ** 2)
     return fac