diff --git a/py4DSTEM/process/polar/polar_analysis.py b/py4DSTEM/process/polar/polar_analysis.py
index 0c2454289..829bbde59 100644
--- a/py4DSTEM/process/polar/polar_analysis.py
+++ b/py4DSTEM/process/polar/polar_analysis.py
@@ -11,15 +11,15 @@
 
 def calculate_radial_statistics(
     self,
-    median_local = False,
-    median_global = False,
-    plot_results_mean = False,
-    plot_results_var = False,
-    figsize = (8,4),
-    returnval = False,
-    returnfig = False,
-    progress_bar = True,
-    ):
+    median_local=False,
+    median_global=False,
+    plot_results_mean=False,
+    plot_results_var=False,
+    figsize=(8, 4),
+    returnval=False,
+    returnfig=False,
+    progress_bar=True,
+):
     """
     Calculate fluctuation electron microscopy (FEM) statistics, including radial mean,
     variance, and normalized variance. This function uses the original FEM definitions,
@@ -49,16 +49,20 @@ def calculate_radial_statistics(
     self.scattering_vector_units = self.calibration.get_Q_pixel_units()
 
     # init radial data arrays
-    self.radial_all = np.zeros((
-        self._datacube.shape[0],
-        self._datacube.shape[1],
-        self.polar_shape[1],
-    ))
-    self.radial_all_std = np.zeros((
-        self._datacube.shape[0],
-        self._datacube.shape[1],
-        self.polar_shape[1],
-    ))
+    self.radial_all = np.zeros(
+        (
+            self._datacube.shape[0],
+            self._datacube.shape[1],
+            self.polar_shape[1],
+        )
+    )
+    self.radial_all_std = np.zeros(
+        (
+            self._datacube.shape[0],
+            self._datacube.shape[1],
+            self.polar_shape[1],
+        )
+    )
 
     # Compute the radial mean for each probe position
     for rx, ry in tqdmnd(
@@ -66,28 +70,26 @@ def calculate_radial_statistics(
         self._datacube.shape[1],
         desc="Radial statistics",
         unit=" probe positions",
-        disable=not progress_bar):
-        
-        self.radial_all[rx,ry] = np.mean(
-            self.data[rx,ry],
-            axis=0)
-        self.radial_all_std[rx,ry] = np.sqrt(np.mean(
-            (self.data[rx,ry] - self.radial_all[rx,ry][None])**2, 
-            axis=0))
-
-    self.radial_mean = np.mean(self.radial_all, axis=(0,1))
+        disable=not progress_bar,
+    ):
+        self.radial_all[rx, ry] = np.mean(self.data[rx, ry], axis=0)
+        self.radial_all_std[rx, ry] = np.sqrt(
+            np.mean((self.data[rx, ry] - self.radial_all[rx, ry][None]) ** 2, axis=0)
+        )
+
+    self.radial_mean = np.mean(self.radial_all, axis=(0, 1))
     self.radial_var = np.mean(
-        (self.radial_all - self.radial_mean[None,None])**2,
-        axis=(0,1))
+        (self.radial_all - self.radial_mean[None, None]) ** 2, axis=(0, 1)
+    )
 
-    self.radial_var_norm = self.radial_var 
+    self.radial_var_norm = self.radial_var
     sub = self.radial_mean > 0.0
-    self.radial_var_norm[sub] /= self.radial_mean[sub]**2
+    self.radial_var_norm[sub] /= self.radial_mean[sub] ** 2
 
     # plot results
     if plot_results_mean:
         if returnfig:
-            fig,ax = plot_radial_mean(
+            fig, ax = plot_radial_mean(
                 self,
                 figsize=figsize,
                 returnfig=True,
@@ -95,20 +97,20 @@ def calculate_radial_statistics(
         else:
             plot_radial_mean(
                 self,
-                figsize = figsize,
-                )
+                figsize=figsize,
+            )
     elif plot_results_var:
         if returnfig:
-            fig,ax = plot_radial_var_norm(
+            fig, ax = plot_radial_var_norm(
                 self,
-                figsize = figsize,
-                returnfig = True,
-                )
+                figsize=figsize,
+                returnfig=True,
+            )
         else:
             plot_radial_var_norm(
                 self,
-                figsize = figsize,
-                )
+                figsize=figsize,
+            )
 
     # Return values
     if returnval:
@@ -125,31 +127,31 @@ def calculate_radial_statistics(
 
 def plot_radial_mean(
     self,
-    log_x = False,
-    log_y = False,
-    figsize = (8,4),
-    returnfig = False,
-    ):
+    log_x=False,
+    log_y=False,
+    figsize=(8, 4),
+    returnfig=False,
+):
     """
     Plot radial mean
     """
-    fig,ax = plt.subplots(figsize=figsize)
+    fig, ax = plt.subplots(figsize=figsize)
     ax.plot(
         self.scattering_vector,
         self.radial_mean,
-        )
+    )
 
     if log_x:
-        ax.set_xscale('log')
+        ax.set_xscale("log")
     if log_y:
-        ax.set_yscale('log')
+        ax.set_yscale("log")
 
-    ax.set_xlabel('Scattering Vector (' + self.scattering_vector_units + ')')
-    ax.set_ylabel('Radial Mean')
+    ax.set_xlabel("Scattering Vector (" + self.scattering_vector_units + ")")
+    ax.set_ylabel("Radial Mean")
     if log_x and self.scattering_vector[0] == 0.0:
-        ax.set_xlim((self.scattering_vector[1],self.scattering_vector[-1]))
+        ax.set_xlim((self.scattering_vector[1], self.scattering_vector[-1]))
     else:
-        ax.set_xlim((self.scattering_vector[0],self.scattering_vector[-1]))
+        ax.set_xlim((self.scattering_vector[0], self.scattering_vector[-1]))
 
     if returnfig:
         return fig, ax
@@ -157,9 +159,9 @@ def plot_radial_mean(
 
 def plot_radial_var_norm(
     self,
-    figsize = (8,4),
-    returnfig = False,
-    ):
+    figsize=(8, 4),
+    returnfig=False,
+):
     """
     Plotting function for the global FEM.
     """
@@ -179,27 +181,27 @@ def plot_radial_var_norm(
 
 def calculate_pair_dist_function(
     self,
-    k_min = 0.05,
-    k_max = None,
-    k_width = 0.25,
-    k_lowpass = None,
-    k_highpass = None,
+    k_min=0.05,
+    k_max=None,
+    k_width=0.25,
+    k_lowpass=None,
+    k_highpass=None,
     # k_pad_max = 10.0,
-    r_min = 0.0,
-    r_max = 20.0,
-    r_step = 0.02,
-    damp_origin_fluctuations = False,
+    r_min=0.0,
+    r_max=20.0,
+    r_step=0.02,
+    damp_origin_fluctuations=False,
     # poly_background_order = 2,
     # iterative_pdf_refine = True,
     # num_iter = 10,
-    dens = None,
-    plot_fits = False,
-    plot_sf_estimate = False,
-    plot_reduced_pdf = True,
-    plot_pdf = False,
-    figsize = (8,4),
-    maxfev = None,
-    ):
+    dens=None,
+    plot_fits=False,
+    plot_sf_estimate=False,
+    plot_reduced_pdf=True,
+    plot_pdf=False,
+    figsize=(8, 4),
+    maxfev=None,
+):
     """
     Calculate the pair distribution function (PDF).
 
@@ -218,7 +220,7 @@ def calculate_pair_dist_function(
     int0 = np.median(self.radial_mean) / int_mean - const_bg
     sigma0 = np.mean(k)
     coefs = [const_bg, int0, sigma0, int0, sigma0]
-    lb = [0,0,0,0,0]
+    lb = [0, 0, 0, 0, 0]
     ub = [np.inf, np.inf, np.inf, np.inf, np.inf]
     # Weight the fit towards high k values
     noise_est = k[-1] - k + dk
@@ -226,31 +228,30 @@ def calculate_pair_dist_function(
     # Estimate the mean atomic form factor + background
     if maxfev is None:
         coefs = curve_fit(
-            scattering_model, 
-            k2[sub_fit], 
-            Ik[sub_fit] / int_mean, 
-            sigma = noise_est[sub_fit],
+            scattering_model,
+            k2[sub_fit],
+            Ik[sub_fit] / int_mean,
+            sigma=noise_est[sub_fit],
             p0=coefs,
-            xtol = 1e-8,
-            bounds = (lb,ub),
+            xtol=1e-8,
+            bounds=(lb, ub),
         )[0]
     else:
         coefs = curve_fit(
-            scattering_model, 
-            k2[sub_fit], 
-            Ik[sub_fit] / int_mean, 
-            sigma = noise_est[sub_fit],
+            scattering_model,
+            k2[sub_fit],
+            Ik[sub_fit] / int_mean,
+            sigma=noise_est[sub_fit],
             p0=coefs,
-            xtol = 1e-8,
-            bounds = (lb,ub),
-            maxfev = maxfev,
+            xtol=1e-8,
+            bounds=(lb, ub),
+            maxfev=maxfev,
         )[0]
 
     coefs[0] *= int_mean
     coefs[1] *= int_mean
     coefs[3] *= int_mean
 
-
     # Calculate the mean atomic form factor wthout any background
     coefs_fk = (0.0, coefs[1], coefs[2], coefs[3], coefs[4])
     fk = scattering_model(k2, coefs_fk)
@@ -263,40 +264,40 @@ def calculate_pair_dist_function(
     #     (k - k_min) / k_width,
     #     (k_max - k) / k_width,
     #     ),0,1)
-    mask = np.clip(np.minimum(
-        (k - 0.0) / k_width,
-        (k_max - k) / k_width,
-        ),0,1)
-    mask = np.sin(mask*(np.pi/2))
+    mask = np.clip(
+        np.minimum(
+            (k - 0.0) / k_width,
+            (k_max - k) / k_width,
+        ),
+        0,
+        1,
+    )
+    mask = np.sin(mask * (np.pi / 2))
 
     # Estimate the reduced structure factor S(k)
     Sk = (Ik - bg) * k / fk
 
     # Masking edges of S(k)
     mask_sum = np.sum(mask)
-    Sk = (Sk - np.sum(Sk*mask)/mask_sum) * mask
+    Sk = (Sk - np.sum(Sk * mask) / mask_sum) * mask
 
     # Filtering of S(k)
     if k_lowpass is not None and k_lowpass > 0.0:
-        Sk = gaussian_filter(
-            Sk, 
-            sigma=k_lowpass / dk,
-            mode = 'nearest')
+        Sk = gaussian_filter(Sk, sigma=k_lowpass / dk, mode="nearest")
     if k_highpass is not None:
-        Sk_lowpass = gaussian_filter(
-            Sk, 
-            sigma=k_highpass / dk,
-            mode = 'nearest')
+        Sk_lowpass = gaussian_filter(Sk, sigma=k_highpass / dk, mode="nearest")
         Sk -= Sk_lowpass
 
     # Calculate the real space PDF
     r = np.arange(r_min, r_max, r_step)
-    ra,ka = np.meshgrid(r,k)
-    pdf_reduced = (2/np.pi)*dk*np.sum(
-        np.sin(
-            2*np.pi*ra*ka
-        ) * Sk[:,None],
-        axis=0,
+    ra, ka = np.meshgrid(r, k)
+    pdf_reduced = (
+        (2 / np.pi)
+        * dk
+        * np.sum(
+            np.sin(2 * np.pi * ra * ka) * Sk[:, None],
+            axis=0,
+        )
     )
 
     # Damp the unphysical fluctuations at the PDF origin
@@ -304,7 +305,7 @@ def calculate_pair_dist_function(
         ind_max = np.argmax(pdf_reduced)
         r_ind_max = r[ind_max]
         r_mask = np.minimum(r / r_ind_max, 1.0)
-        r_mask = np.sin(r_mask*np.pi/2)**2
+        r_mask = np.sin(r_mask * np.pi / 2) ** 2
         pdf_reduced *= r_mask
 
     # Store results
@@ -314,8 +315,8 @@ def calculate_pair_dist_function(
     # if density is provided, we can estimate the full PDF
     if dens is not None:
         pdf = pdf_reduced.copy()
-        pdf[1:] /= (4*np.pi*dens*r[1:]*(r[1]-r[0]))
-        pdf *= (2/np.pi)
+        pdf[1:] /= 4 * np.pi * dens * r[1:] * (r[1] - r[0])
+        pdf *= 2 / np.pi
         pdf += 1
 
         if damp_origin_fluctuations:
@@ -323,68 +324,70 @@ def calculate_pair_dist_function(
 
         pdf = np.maximum(pdf, 0.0)
 
-
-
     # Plots
     if plot_fits:
-        fig,ax = plt.subplots(figsize=figsize)
+        fig, ax = plt.subplots(figsize=figsize)
         ax.plot(
             self.scattering_vector,
             self.radial_mean,
-            color = 'k',
-            )
+            color="k",
+        )
         ax.plot(
             k,
             bg,
-            color = 'r',
-            )
-        ax.set_xlabel('Scattering Vector (' + self.scattering_vector_units + ')')
-        ax.set_ylabel('Radial Mean')
-        ax.set_xlim((self.scattering_vector[0],self.scattering_vector[-1]))
+            color="r",
+        )
+        ax.set_xlabel("Scattering Vector (" + self.scattering_vector_units + ")")
+        ax.set_ylabel("Radial Mean")
+        ax.set_xlim((self.scattering_vector[0], self.scattering_vector[-1]))
         # ax.set_ylim((0,2e-5))
-        ax.set_xlabel('Scattering Vector [A^-1]')
-        ax.set_ylabel('I(k) and Fit Estimates')
+        ax.set_xlabel("Scattering Vector [A^-1]")
+        ax.set_ylabel("I(k) and Fit Estimates")
 
-        ax.set_ylim((np.min(self.radial_mean[self.radial_mean>0])*0.8,
-            np.max(self.radial_mean*mask)*1.25))
-        ax.set_yscale('log')
+        ax.set_ylim(
+            (
+                np.min(self.radial_mean[self.radial_mean > 0]) * 0.8,
+                np.max(self.radial_mean * mask) * 1.25,
+            )
+        )
+        ax.set_yscale("log")
 
     if plot_sf_estimate:
-        fig,ax = plt.subplots(figsize=figsize)
+        fig, ax = plt.subplots(figsize=figsize)
         ax.plot(
             k,
             Sk,
-            color = 'r',
+            color="r",
+        )
+        yr = (np.min(Sk), np.max(Sk))
+        ax.set_ylim(
+            (
+                yr[0] - 0.05 * (yr[1] - yr[0]),
+                yr[1] + 0.05 * (yr[1] - yr[0]),
             )
-        yr = (np.min(Sk),np.max(Sk))
-        ax.set_ylim((
-            yr[0]-0.05*(yr[1]-yr[0]),
-            yr[1]+0.05*(yr[1]-yr[0]),
-            ))
-        ax.set_xlabel('Scattering Vector [A^-1]')
-        ax.set_ylabel('Reduced Structure Factor')
+        )
+        ax.set_xlabel("Scattering Vector [A^-1]")
+        ax.set_ylabel("Reduced Structure Factor")
 
     if plot_reduced_pdf:
-        fig,ax = plt.subplots(figsize=figsize)
+        fig, ax = plt.subplots(figsize=figsize)
         ax.plot(
             r,
             pdf_reduced,
-            color = 'r',
-            )
-        ax.set_xlabel('Radius [A]')
-        ax.set_ylabel('Reduced Pair Distribution Function')
+            color="r",
+        )
+        ax.set_xlabel("Radius [A]")
+        ax.set_ylabel("Reduced Pair Distribution Function")
 
     if plot_pdf:
-        fig,ax = plt.subplots(figsize=figsize)
+        fig, ax = plt.subplots(figsize=figsize)
         ax.plot(
             r,
             pdf,
-            color = 'r',
-            )
-        ax.set_xlabel('Radius [A]')
-        ax.set_ylabel('Pair Distribution Function')
-
-
+            color="r",
+        )
+        ax.set_xlabel("Radius [A]")
+        ax.set_ylabel("Pair Distribution Function")
 
     # functions for inverting from reduced PDF back to S(k)
 
@@ -403,7 +406,6 @@ def calculate_pair_dist_function(
     # )
 
 
-
 def calculate_FEM_local(
     self,
     figsize=(8, 6),
@@ -428,13 +430,10 @@ def calculate_FEM_local(
 
 
     """
-    
 
     pass
 
 
-
-
 def scattering_model(k2, *coefs):
     coefs = np.squeeze(np.array(coefs))
 
@@ -444,14 +443,14 @@ def scattering_model(k2, *coefs):
     int1 = coefs[3]
     sigma1 = coefs[4]
 
-    int_model = const_bg + \
-        int0*np.exp(k2/(-2*sigma0**2)) + \
-        int1*np.exp(k2**2/(-2*sigma1**4))
-
-        # (int1*sigma1)/(k2 + sigma1**2)
-        # int1*np.exp(k2/(-2*sigma1**2))
-        # int1*np.exp(k2/(-2*sigma1**2))
+    int_model = (
+        const_bg
+        + int0 * np.exp(k2 / (-2 * sigma0**2))
+        + int1 * np.exp(k2**2 / (-2 * sigma1**4))
+    )
 
+    # (int1*sigma1)/(k2 + sigma1**2)
+    # int1*np.exp(k2/(-2*sigma1**2))
+    # int1*np.exp(k2/(-2*sigma1**2))
 
     return int_model
-