diff --git a/package/lib/cross_correlation.py b/package/lib/cross_correlation.py
index 5613c15..159073d 100755
--- a/package/lib/cross_correlation.py
+++ b/package/lib/cross_correlation.py
@@ -48,20 +48,16 @@ def _upsampled_dft(data, upsampled_region_size, upsample_factor=1, axis_offsets=
     """
     # if people pass in an integer, expand it to a list of equal-sized sections
     if not hasattr(upsampled_region_size, "__iter__"):
-        upsampled_region_size = [
-            upsampled_region_size,
-        ] * data.ndim
+        upsampled_region_size = [upsampled_region_size] * data.ndim
     else:
         if len(upsampled_region_size) != data.ndim:
-            raise ValueError("shape of upsampled region sizes must be equal " "to input data's number of dimensions.")
+            raise ValueError("shape of upsampled region sizes must be equal to input data's number of dimensions.")
 
     if axis_offsets is None:
-        axis_offsets = [
-            0,
-        ] * data.ndim
+        axis_offsets = [0] * data.ndim
     else:
         if len(axis_offsets) != data.ndim:
-            raise ValueError("number of axis offsets must be equal to input " "data's number of dimensions.")
+            raise ValueError("number of axis offsets must be equal to input data's number of dimensions.")
 
     im2pi = 1j * 2 * np.pi
 
diff --git a/package/lib/plots.py b/package/lib/plots.py
index 40ec567..c5a92d2 100755
--- a/package/lib/plots.py
+++ b/package/lib/plots.py
@@ -54,10 +54,7 @@ from matplotlib.colors import LogNorm
 from matplotlib.patches import Circle, FancyArrowPatch, Rectangle
 from matplotlib.path import Path
 from matplotlib.widgets import Button, LassoSelector, RectangleSelector, Slider, TextBox
-from mpl_toolkits.axes_grid1.anchored_artists import (
-    AnchoredDirectionArrows,
-    AnchoredSizeBar,
-)
+from mpl_toolkits.axes_grid1.anchored_artists import AnchoredDirectionArrows, AnchoredSizeBar
 from scipy.ndimage import zoom as sc_zoom
 
 try:
@@ -66,15 +63,7 @@ except ImportError:
     from utils import PCconf, princ_angle, rot2D, sci_not
 
 
-def plot_obs(
-    data_array,
-    headers,
-    rectangle=None,
-    shifts=None,
-    savename=None,
-    plots_folder="",
-    **kwargs,
-):
+def plot_obs(data_array, headers, rectangle=None, shifts=None, savename=None, plots_folder="", **kwargs):
     """
     Plots raw observation imagery with some information on the instrument and
     filters.
@@ -105,14 +94,7 @@ def plot_obs(
     plt.rcParams.update({"font.size": 10})
     nb_obs = np.max([np.sum([head["filtnam1"] == curr_pol for head in headers]) for curr_pol in ["POL0", "POL60", "POL120"]])
     shape = np.array((3, nb_obs))
-    fig, ax = plt.subplots(
-        shape[0],
-        shape[1],
-        figsize=(3 * shape[1], 3 * shape[0]),
-        layout="compressed",
-        sharex=True,
-        sharey=True,
-    )
+    fig, ax = plt.subplots(shape[0], shape[1], figsize=(3 * shape[1], 3 * shape[0]), layout="compressed", sharex=True, sharey=True)
     used = np.zeros(shape, dtype=bool)
     r_pol = dict(pol0=0, pol60=1, pol120=2)
     c_pol = dict(pol0=0, pol60=0, pol120=0)
@@ -135,14 +117,8 @@ def plot_obs(
             vmin, vmax = kwargs["vmin"], kwargs["vmax"]
             del kwargs["vmin"], kwargs["vmax"]
         else:
-            vmin, vmax = (
-                convert * data[data > 0.0].min() / 10.0,
-                convert * data[data > 0.0].max(),
-            )
-        for key, value in [
-            ["cmap", [["cmap", "inferno"]]],
-            ["norm", [["norm", LogNorm(vmin, vmax)]]],
-        ]:
+            vmin, vmax = (convert * data[data > 0.0].min() / 10.0, convert * data[data > 0.0].max())
+        for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]:
             try:
                 _ = kwargs[key]
             except KeyError:
@@ -154,64 +130,21 @@ def plot_obs(
         if shifts is not None:
             x, y = np.array(data.shape[::-1]) / 2.0 - shifts[i]
             dx, dy = shifts[i]
-            ax_curr.arrow(
-                x,
-                y,
-                dx,
-                dy,
-                length_includes_head=True,
-                width=0.1,
-                head_width=0.3,
-                color="g",
-            )
+            ax_curr.arrow(x, y, dx, dy, length_includes_head=True, width=0.1, head_width=0.3, color="g")
             ax_curr.plot([x, x], [0, data.shape[0] - 1], "--", lw=2, color="g", alpha=0.85)
             ax_curr.plot([0, data.shape[1] - 1], [y, y], "--", lw=2, color="g", alpha=0.85)
         if rectangle is not None:
             x, y, width, height, angle, color = rectangle[i]
             ax_curr.add_patch(Rectangle((x, y), width, height, angle=angle, edgecolor=color, fill=False))
         # position of centroid
-        ax_curr.plot(
-            [data.shape[1] / 2, data.shape[1] / 2],
-            [0, data.shape[0] - 1],
-            "--",
-            lw=2,
-            color="b",
-            alpha=0.85,
-        )
-        ax_curr.plot(
-            [0, data.shape[1] - 1],
-            [data.shape[0] / 2, data.shape[0] / 2],
-            "--",
-            lw=2,
-            color="b",
-            alpha=0.85,
-        )
+        ax_curr.plot([data.shape[1] / 2, data.shape[1] / 2], [0, data.shape[0] - 1], "--", lw=2, color="b", alpha=0.85)
+        ax_curr.plot([0, data.shape[1] - 1], [data.shape[0] / 2, data.shape[0] / 2], "--", lw=2, color="b", alpha=0.85)
         ax_curr.annotate(
-            instr + ":" + rootname,
-            color="white",
-            fontsize=10,
-            xy=(0.01, 1.00),
-            xycoords="axes fraction",
-            verticalalignment="top",
-            horizontalalignment="left",
+            instr + ":" + rootname, color="white", fontsize=10, xy=(0.01, 1.00), xycoords="axes fraction", verticalalignment="top", horizontalalignment="left"
         )
+        ax_curr.annotate(filt, color="white", fontsize=15, xy=(0.01, 0.01), xycoords="axes fraction", verticalalignment="bottom", horizontalalignment="left")
         ax_curr.annotate(
-            filt,
-            color="white",
-            fontsize=15,
-            xy=(0.01, 0.01),
-            xycoords="axes fraction",
-            verticalalignment="bottom",
-            horizontalalignment="left",
-        )
-        ax_curr.annotate(
-            exptime,
-            color="white",
-            fontsize=10,
-            xy=(1.00, 0.01),
-            xycoords="axes fraction",
-            verticalalignment="bottom",
-            horizontalalignment="right",
+            exptime, color="white", fontsize=10, xy=(1.00, 0.01), xycoords="axes fraction", verticalalignment="bottom", horizontalalignment="right"
         )
     unused = np.logical_not(used)
     ii, jj = np.indices(shape)
@@ -219,26 +152,13 @@ def plot_obs(
         fig.delaxes(ax[i][j])
 
     # fig.subplots_adjust(hspace=0.01, wspace=0.01, right=1.02)
-    fig.colorbar(
-        im,
-        ax=ax,
-        location="right",
-        shrink=0.75,
-        aspect=50,
-        pad=0.025,
-        label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]",
-    )
+    fig.colorbar(im, ax=ax, location="right", shrink=0.75, aspect=50, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
 
     if savename is not None:
         # fig.suptitle(savename)
         if savename[-4:] not in [".png", ".jpg", ".pdf"]:
             savename += ".pdf"
-        fig.savefig(
-            path_join(plots_folder, savename),
-            bbox_inches="tight",
-            dpi=150,
-            facecolor="None",
-        )
+        fig.savefig(path_join(plots_folder, savename), bbox_inches="tight", dpi=150, facecolor="None")
     plt.show()
     return 0
 
@@ -297,12 +217,7 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""):
             savename += "_IQU.pdf"
         else:
             savename = savename[:-4] + "_IQU" + savename[-4:]
-        fig.savefig(
-            path_join(plots_folder, savename),
-            bbox_inches="tight",
-            dpi=150,
-            facecolor="None",
-        )
+        fig.savefig(path_join(plots_folder, savename), bbox_inches="tight", dpi=150, facecolor="None")
     return 0
 
 
@@ -389,10 +304,7 @@ def polarization_map(
 
     # Compute confidence level map
     QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan)
-    for nflux, sflux in zip(
-        [QN, UN, QN_ERR, UN_ERR],
-        [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])],
-    ):
+    for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]):
         nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0]
     confP = PCconf(QN, UN, QN_ERR, UN_ERR)
 
@@ -461,11 +373,7 @@ def polarization_map(
     ax.set_ylabel("Declination (J2000)", labelpad=-1)
 
     vmin, vmax = 0.0, stkI.max() * convert_flux
-    for key, value in [
-        ["cmap", [["cmap", "inferno"]]],
-        ["width", [["width", 0.4]]],
-        ["linewidth", [["linewidth", 0.6]]],
-    ]:
+    for key, value in [["cmap", [["cmap", "inferno"]]], ["width", [["width", 0.4]]], ["linewidth", [["linewidth", 0.6]]]]:
         try:
             _ = kwargs[key]
         except KeyError:
@@ -623,18 +531,7 @@ def polarization_map(
 
     plt.rcParams.update({"font.size": 11})
     px_size = wcs.wcs.get_cdelt()[0] * 3600.0
-    px_sc = AnchoredSizeBar(
-        ax.transData,
-        1.0 / px_size,
-        "1 arcsec",
-        3,
-        pad=0.25,
-        sep=5,
-        borderpad=0.25,
-        frameon=False,
-        size_vertical=0.005,
-        color=font_color,
-    )
+    px_sc = AnchoredSizeBar(ax.transData, 1.0 / px_size, "1 arcsec", 3, pad=0.25, sep=5, borderpad=0.25, frameon=False, size_vertical=0.005, color=font_color)
     north_dir = AnchoredDirectionArrows(
         ax.transAxes,
         "E",
@@ -659,10 +556,7 @@ def polarization_map(
             step_vec = 1
             scale_vec = 2.0
         X, Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
-        U, V = (
-            poldata * np.cos(np.pi / 2.0 + pangdata * np.pi / 180.0),
-            poldata * np.sin(np.pi / 2.0 + pangdata * np.pi / 180.0),
-        )
+        U, V = (poldata * np.cos(np.pi / 2.0 + pangdata * np.pi / 180.0), poldata * np.sin(np.pi / 2.0 + pangdata * np.pi / 180.0))
         ax.quiver(
             X[::step_vec, ::step_vec],
             Y[::step_vec, ::step_vec],
@@ -682,16 +576,7 @@ def polarization_map(
             edgecolor="k",
         )
         pol_sc = AnchoredSizeBar(
-            ax.transData,
-            scale_vec,
-            r"$P$= 100 %",
-            4,
-            pad=0.25,
-            sep=5,
-            borderpad=0.25,
-            frameon=False,
-            size_vertical=0.005,
-            color=font_color,
+            ax.transData, scale_vec, r"$P$= 100 %", 4, pad=0.25, sep=5, borderpad=0.25, frameon=False, size_vertical=0.005, color=font_color
         )
 
         ax.add_artist(pol_sc)
@@ -700,8 +585,7 @@ def polarization_map(
 
         ax.annotate(
             r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(
-                pivot_wav,
-                sci_not(I_diluted * convert_flux, I_diluted_err * convert_flux, 2),
+                pivot_wav, sci_not(I_diluted * convert_flux, I_diluted_err * convert_flux, 2)
             )
             + "\n"
             + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_diluted * 100.0, P_diluted_err * 100.0)
@@ -720,8 +604,7 @@ def polarization_map(
             ax.add_artist(north_dir)
         ax.annotate(
             r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(
-                pivot_wav,
-                sci_not(I_diluted * convert_flux, I_diluted_err * convert_flux, 2),
+                pivot_wav, sci_not(I_diluted * convert_flux, I_diluted_err * convert_flux, 2)
             ),
             color="white",
             xy=(0.01, 1.00),
@@ -740,12 +623,7 @@ def polarization_map(
     if savename is not None:
         if savename[-4:] not in [".png", ".jpg", ".pdf"]:
             savename += ".pdf"
-        fig.savefig(
-            path_join(plots_folder, savename),
-            bbox_inches="tight",
-            dpi=300,
-            facecolor="None",
-        )
+        fig.savefig(path_join(plots_folder, savename), bbox_inches="tight", dpi=300, facecolor="None")
 
     return fig, ax
 
@@ -781,26 +659,14 @@ class align_maps(object):
             self.other_data = self.other_data[0]
 
         self.map_convert, self.map_unit = (
-            (
-                float(self.map_header["photflam"]),
-                r"$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$",
-            )
+            (float(self.map_header["photflam"]), r"$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$")
             if "PHOTFLAM" in list(self.map_header.keys())
-            else (
-                1.0,
-                self.map_header["bunit"] if "BUNIT" in list(self.map_header.keys()) else "Arbitray Units",
-            )
+            else (1.0, self.map_header["bunit"] if "BUNIT" in list(self.map_header.keys()) else "Arbitray Units")
         )
         self.other_convert, self.other_unit = (
-            (
-                float(self.other_header["photflam"]),
-                r"$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$",
-            )
+            (float(self.other_header["photflam"]), r"$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$")
             if "PHOTFLAM" in list(self.other_header.keys())
-            else (
-                1.0,
-                self.other_header["bunit"] if "BUNIT" in list(self.other_header.keys()) else "Arbitray Units",
-            )
+            else (1.0, self.other_header["bunit"] if "BUNIT" in list(self.other_header.keys()) else "Arbitray Units")
         )
         self.map_observer = (
             "/".join([self.map_header["telescop"], self.map_header["instrume"]]) if "INSTRUME" in list(self.map_header.keys()) else self.map_header["telescop"]
@@ -819,14 +685,8 @@ class align_maps(object):
 
         # Plot the UV map
         other_kwargs = deepcopy(kwargs)
-        vmin, vmax = (
-            self.map_data[self.map_data > 0.0].max() / 1e3 * self.map_convert,
-            self.map_data[self.map_data > 0.0].max() * self.map_convert,
-        )
-        for key, value in [
-            ["cmap", [["cmap", "inferno"]]],
-            ["norm", [["norm", LogNorm(vmin, vmax)]]],
-        ]:
+        vmin, vmax = (self.map_data[self.map_data > 0.0].max() / 1e3 * self.map_convert, self.map_data[self.map_data > 0.0].max() * self.map_convert)
+        for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]:
             try:
                 _ = kwargs[key]
             except KeyError:
@@ -918,10 +778,7 @@ class align_maps(object):
             self.other_data[self.other_data > 0.0].max() / 1e3 * self.other_convert,
             self.other_data[self.other_data > 0.0].max() * self.other_convert,
         )
-        for key, value in [
-            ["cmap", [["cmap", "inferno"]]],
-            ["norm", [["norm", LogNorm(vmin, vmax)]]],
-        ]:
+        for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]:
             try:
                 _ = other_kwargs[key]
             except KeyError:
@@ -1025,10 +882,7 @@ class align_maps(object):
         else:
             self.map_wcs.wcs.crpix = np.array(self.cr_map.get_data()) + (1.0, 1.0)
         if np.array(self.cr_other.get_data()).shape == (2, 1):
-            self.other_wcs.wcs.crpix = np.array(self.cr_other.get_data())[:, 0] + (
-                1.0,
-                1.0,
-            )
+            self.other_wcs.wcs.crpix = np.array(self.cr_other.get_data())[:, 0] + (1.0, 1.0)
         else:
             self.other_wcs.wcs.crpix = np.array(self.cr_other.get_data()) + (1.0, 1.0)
         self.map_wcs.wcs.crval = np.array(self.map_wcs.pixel_to_world_values(*self.map_wcs.wcs.crpix))
@@ -1080,15 +934,7 @@ class overplot_radio(align_maps):
     Inherit from class align_maps in order to get the same WCS on both maps.
     """
 
-    def overplot(
-        self,
-        levels=None,
-        P_cut=0.99,
-        SNRi_cut=1.0,
-        scale_vec=2,
-        savename=None,
-        **kwargs,
-    ):
+    def overplot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, scale_vec=2, savename=None, **kwargs):
         self.Stokes_UV = self.map
         self.wcs_UV = self.map_wcs
         # Get Data
@@ -1102,10 +948,7 @@ class overplot_radio(align_maps):
         pang = self.Stokes_UV["POL_ANG"].data
         # Compute confidence level map
         QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan)
-        for nflux, sflux in zip(
-            [QN, UN, QN_ERR, UN_ERR],
-            [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])],
-        ):
+        for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]):
             nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0]
         confP = PCconf(QN, UN, QN_ERR, UN_ERR)
 
@@ -1137,14 +980,8 @@ class overplot_radio(align_maps):
         self.fig_overplot.subplots_adjust(hspace=0, wspace=0, bottom=0.1, left=0.1, top=0.8, right=1)
 
         # Display UV intensity map with polarization vectors
-        vmin, vmax = (
-            stkI[np.isfinite(stkI)].max() / 1e3 * self.map_convert,
-            stkI[np.isfinite(stkI)].max() * self.map_convert,
-        )
-        for key, value in [
-            ["cmap", [["cmap", "inferno"]]],
-            ["norm", [["norm", LogNorm(vmin, vmax)]]],
-        ]:
+        vmin, vmax = (stkI[np.isfinite(stkI)].max() / 1e3 * self.map_convert, stkI[np.isfinite(stkI)].max() * self.map_convert)
+        for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]:
             try:
                 _ = kwargs[key]
             except KeyError:
@@ -1179,19 +1016,9 @@ class overplot_radio(align_maps):
         else:
             self.ax_overplot.set_facecolor("white")
             font_color = "black"
-        self.im = self.ax_overplot.imshow(
-            stkI * self.map_convert,
-            aspect="equal",
-            label="{0:s} observation".format(self.map_observer),
-            **kwargs,
-        )
+        self.im = self.ax_overplot.imshow(stkI * self.map_convert, aspect="equal", label="{0:s} observation".format(self.map_observer), **kwargs)
         self.cbar = self.fig_overplot.colorbar(
-            self.im,
-            ax=self.ax_overplot,
-            aspect=50,
-            shrink=0.75,
-            pad=0.025,
-            label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit),
+            self.im, ax=self.ax_overplot, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit)
         )
 
         # Display full size polarization vectors
@@ -1202,10 +1029,7 @@ class overplot_radio(align_maps):
             self.scale_vec = scale_vec
         step_vec = 1
         self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
-        self.U, self.V = (
-            pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0),
-            pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0),
-        )
+        self.U, self.V = (pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0))
         self.Q = self.ax_overplot.quiver(
             self.X[::step_vec, ::step_vec],
             self.Y[::step_vec, ::step_vec],
@@ -1242,11 +1066,7 @@ class overplot_radio(align_maps):
         self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1)
         self.fig_overplot.suptitle(
             "{0:s} polarization map of {1:s} overplotted with {2:s} {3:.2f}GHz map in {4:s}.".format(
-                self.map_observer,
-                obj,
-                self.other_observer,
-                other_freq * 1e-9,
-                self.other_unit,
+                self.map_observer, obj, self.other_observer, other_freq * 1e-9, self.other_unit
             ),
             wrap=True,
         )
@@ -1300,9 +1120,7 @@ class overplot_radio(align_maps):
 
         (self.cr_map,) = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix - (1.0, 1.0)), "r+")
         (self.cr_other,) = self.ax_overplot.plot(
-            *(self.other_wcs.celestial.wcs.crpix - (1.0, 1.0)),
-            "g+",
-            transform=self.ax_overplot.get_transform(self.other_wcs),
+            *(self.other_wcs.celestial.wcs.crpix - (1.0, 1.0)), "g+", transform=self.ax_overplot.get_transform(self.other_wcs)
         )
 
         handles, labels = self.ax_overplot.get_legend_handles_labels()
@@ -1312,12 +1130,7 @@ class overplot_radio(align_maps):
         labels.append("{0:s} contour".format(self.other_observer))
         handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.get_edgecolor()[0]))
         self.legend = self.ax_overplot.legend(
-            handles=handles,
-            labels=labels,
-            bbox_to_anchor=(0.0, 1.02, 1.0, 0.102),
-            loc="lower left",
-            mode="expand",
-            borderaxespad=0.0,
+            handles=handles, labels=labels, bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), loc="lower left", mode="expand", borderaxespad=0.0
         )
 
         if savename is not None:
@@ -1340,16 +1153,7 @@ class overplot_chandra(align_maps):
     Inherit from class align_maps in order to get the same WCS on both maps.
     """
 
-    def overplot(
-        self,
-        levels=None,
-        P_cut=0.99,
-        SNRi_cut=1.0,
-        scale_vec=2,
-        zoom=1,
-        savename=None,
-        **kwargs,
-    ):
+    def overplot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, scale_vec=2, zoom=1, savename=None, **kwargs):
         self.Stokes_UV = self.map
         self.wcs_UV = self.map_wcs
         # Get Data
@@ -1363,10 +1167,7 @@ class overplot_chandra(align_maps):
         pang = self.Stokes_UV["POL_ANG"].data
         # Compute confidence level map
         QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan)
-        for nflux, sflux in zip(
-            [QN, UN, QN_ERR, UN_ERR],
-            [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])],
-        ):
+        for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]):
             nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0]
         confP = PCconf(QN, UN, QN_ERR, UN_ERR)
 
@@ -1397,14 +1198,8 @@ class overplot_chandra(align_maps):
         self.fig_overplot.subplots_adjust(hspace=0, wspace=0, bottom=0.1, left=0.1, top=0.8, right=1)
 
         # Display UV intensity map with polarization vectors
-        vmin, vmax = (
-            stkI[np.isfinite(stkI)].max() / 1e3 * self.map_convert,
-            stkI[np.isfinite(stkI)].max() * self.map_convert,
-        )
-        for key, value in [
-            ["cmap", [["cmap", "inferno"]]],
-            ["norm", [["norm", LogNorm(vmin, vmax)]]],
-        ]:
+        vmin, vmax = (stkI[np.isfinite(stkI)].max() / 1e3 * self.map_convert, stkI[np.isfinite(stkI)].max() * self.map_convert)
+        for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]:
             try:
                 _ = kwargs[key]
             except KeyError:
@@ -1441,12 +1236,7 @@ class overplot_chandra(align_maps):
             font_color = "black"
         self.im = self.ax_overplot.imshow(stkI * self.map_convert, aspect="equal", **kwargs)
         self.cbar = self.fig_overplot.colorbar(
-            self.im,
-            ax=self.ax_overplot,
-            aspect=50,
-            shrink=0.75,
-            pad=0.025,
-            label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit),
+            self.im, ax=self.ax_overplot, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit)
         )
 
         # Display full size polarization vectors
@@ -1457,10 +1247,7 @@ class overplot_chandra(align_maps):
             self.scale_vec = scale_vec
         step_vec = 1
         self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
-        self.U, self.V = (
-            pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0),
-            pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0),
-        )
+        self.U, self.V = (pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0))
         self.Q = self.ax_overplot.quiver(
             self.X[::step_vec, ::step_vec],
             self.Y[::step_vec, ::step_vec],
@@ -1488,18 +1275,14 @@ class overplot_chandra(align_maps):
         elif zoom != 1:
             levels *= other_data.max() / self.other_data.max()
         other_cont = self.ax_overplot.contour(
-            other_data * self.other_convert,
-            transform=self.ax_overplot.get_transform(other_wcs),
-            levels=levels,
-            colors="grey",
+            other_data * self.other_convert, transform=self.ax_overplot.get_transform(other_wcs), levels=levels, colors="grey"
         )
         self.ax_overplot.clabel(other_cont, inline=True, fontsize=8)
 
         self.ax_overplot.set_xlabel(label="Right Ascension (J2000)")
         self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1)
         self.fig_overplot.suptitle(
-            "{0:s} polarization map of {1:s} overplotted\nwith {2:s} contour in counts.".format(self.map_observer, obj, self.other_observer),
-            wrap=True,
+            "{0:s} polarization map of {1:s} overplotted\nwith {2:s} contour in counts.".format(self.map_observer, obj, self.other_observer), wrap=True
         )
 
         # Display pixel scale and North direction
@@ -1550,11 +1333,7 @@ class overplot_chandra(align_maps):
         self.ax_overplot.add_artist(pol_sc)
 
         (self.cr_map,) = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix - (1.0, 1.0)), "r+")
-        (self.cr_other,) = self.ax_overplot.plot(
-            *(other_wcs.celestial.wcs.crpix - (1.0, 1.0)),
-            "g+",
-            transform=self.ax_overplot.get_transform(other_wcs),
-        )
+        (self.cr_other,) = self.ax_overplot.plot(*(other_wcs.celestial.wcs.crpix - (1.0, 1.0)), "g+", transform=self.ax_overplot.get_transform(other_wcs))
         handles, labels = self.ax_overplot.get_legend_handles_labels()
         handles[np.argmax([li == "{0:s} polarization map".format(self.map_observer) for li in labels])] = FancyArrowPatch(
             (0, 0), (0, 1), arrowstyle="-", fc="w", ec="k", lw=2
@@ -1562,12 +1341,7 @@ class overplot_chandra(align_maps):
         labels.append("{0:s} contour in counts".format(self.other_observer))
         handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.get_edgecolor()[0]))
         self.legend = self.ax_overplot.legend(
-            handles=handles,
-            labels=labels,
-            bbox_to_anchor=(0.0, 1.02, 1.0, 0.102),
-            loc="lower left",
-            mode="expand",
-            borderaxespad=0.0,
+            handles=handles, labels=labels, bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), loc="lower left", mode="expand", borderaxespad=0.0
         )
 
         if savename is not None:
@@ -1580,14 +1354,7 @@ class overplot_chandra(align_maps):
     def plot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, zoom=1, savename=None, **kwargs) -> None:
         while not self.aligned:
             self.align()
-        self.overplot(
-            levels=levels,
-            P_cut=P_cut,
-            SNRi_cut=SNRi_cut,
-            zoom=zoom,
-            savename=savename,
-            **kwargs,
-        )
+        self.overplot(levels=levels, P_cut=P_cut, SNRi_cut=SNRi_cut, zoom=zoom, savename=savename, **kwargs)
         plt.show(block=True)
 
 
@@ -1597,17 +1364,7 @@ class overplot_pol(align_maps):
     Inherit from class align_maps in order to get the same WCS on both maps.
     """
 
-    def overplot(
-        self,
-        levels="Default",
-        P_cut=0.99,
-        SNRi_cut=1.0,
-        step_vec=1,
-        scale_vec=2.0,
-        disptype="i",
-        savename=None,
-        **kwargs,
-    ):
+    def overplot(self, levels="Default", P_cut=0.99, SNRi_cut=1.0, step_vec=1, scale_vec=2.0, disptype="i", savename=None, **kwargs):
         self.Stokes_UV = self.map
         self.wcs_UV = self.map_wcs
         # Get Data
@@ -1621,10 +1378,7 @@ class overplot_pol(align_maps):
         pang = self.Stokes_UV["POL_ANG"].data
         # Compute confidence level map
         QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan)
-        for nflux, sflux in zip(
-            [QN, UN, QN_ERR, UN_ERR],
-            [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])],
-        ):
+        for nflux, sflux in zip([QN, UN, QN_ERR, UN_ERR], [stkQ, stkU, np.sqrt(stk_cov[1, 1]), np.sqrt(stk_cov[2, 2])]):
             nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0]
         confP = PCconf(QN, UN, QN_ERR, UN_ERR)
 
@@ -1648,19 +1402,13 @@ class overplot_pol(align_maps):
         ratiox = max(int(stkI.shape[1] / stkI.shape[0]), 1)
         ratioy = max(int(stkI.shape[0] / stkI.shape[1]), 1)
         self.px_scale = self.wcs_UV.wcs.get_cdelt()[0] / self.other_wcs.wcs.get_cdelt()[0]
-        self.fig_overplot, self.ax_overplot = plt.subplots(
-            figsize=(10 * ratiox, 12 * ratioy),
-            subplot_kw=dict(projection=self.other_wcs),
-        )
+        self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(10 * ratiox, 12 * ratioy), subplot_kw=dict(projection=self.other_wcs))
         self.fig_overplot.subplots_adjust(hspace=0, wspace=0, bottom=0.1, left=0.1, top=0.80, right=1.02)
 
         self.ax_overplot.set_xlabel(label="Right Ascension (J2000)")
         self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1)
         # Display "other" intensity map
-        vmin, vmax = (
-            other_data[other_data > 0.0].max() / 1e3 * self.other_convert,
-            other_data[other_data > 0.0].max() * self.other_convert,
-        )
+        vmin, vmax = (other_data[other_data > 0.0].max() / 1e3 * self.other_convert, other_data[other_data > 0.0].max() * self.other_convert)
         for key, value in [
             ["cmap", [["cmap", "inferno"]]],
             ["norm", [["vmin", vmin], ["vmax", vmax]]],
@@ -1703,11 +1451,7 @@ class overplot_pol(align_maps):
             font_color = "black"
         if "norm" in kwargs.keys():
             self.im = self.ax_overplot.imshow(
-                other_data * self.other_convert,
-                alpha=1.0,
-                label="{0:s} observation".format(self.other_observer),
-                cmap=kwargs["cmap"],
-                norm=kwargs["norm"],
+                other_data * self.other_convert, alpha=1.0, label="{0:s} observation".format(self.other_observer), cmap=kwargs["cmap"], norm=kwargs["norm"]
             )
         else:
             self.im = self.ax_overplot.imshow(
@@ -1719,12 +1463,7 @@ class overplot_pol(align_maps):
                 vmax=kwargs["vmax"],
             )
         self.cbar = self.fig_overplot.colorbar(
-            self.im,
-            ax=self.ax_overplot,
-            aspect=80,
-            shrink=0.75,
-            pad=0.025,
-            label=r"$F_{{\lambda}}$ [{0:s}]".format(self.other_unit),
+            self.im, ax=self.ax_overplot, aspect=80, shrink=0.75, pad=0.025, label=r"$F_{{\lambda}}$ [{0:s}]".format(self.other_unit)
         )
 
         # Display full size polarization vectors
@@ -1737,10 +1476,7 @@ class overplot_pol(align_maps):
             else:
                 self.scale_vec = scale_vec * self.px_scale
             self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
-            self.U, self.V = (
-                pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0),
-                pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0),
-            )
+            self.U, self.V = (pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0))
             self.Q = self.ax_overplot.quiver(
                 self.X[::step_vec, ::step_vec],
                 self.Y[::step_vec, ::step_vec],
@@ -1770,44 +1506,26 @@ class overplot_pol(align_maps):
                 if levels == "Default":
                     levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(pol[stkI > 0.0])
                 cont_stk = self.ax_overplot.contour(
-                    pol * 100.0,
-                    levels=levels * 100.0,
-                    colors="grey",
-                    alpha=0.75,
-                    transform=self.ax_overplot.get_transform(self.wcs_UV),
+                    pol * 100.0, levels=levels * 100.0, colors="grey", alpha=0.75, transform=self.ax_overplot.get_transform(self.wcs_UV)
                 )
             elif disptype.lower() == "pf":
                 disptypestr = "polarized flux"
                 if levels == "Default":
                     levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(stkI[stkI > 0.0] * pol[stkI > 0.0]) * self.map_convert
                 cont_stk = self.ax_overplot.contour(
-                    stkI * pol * self.map_convert,
-                    levels=levels,
-                    colors="grey",
-                    alpha=0.75,
-                    transform=self.ax_overplot.get_transform(self.wcs_UV),
+                    stkI * pol * self.map_convert, levels=levels, colors="grey", alpha=0.75, transform=self.ax_overplot.get_transform(self.wcs_UV)
                 )
             elif disptype.lower() == "snri":
                 disptypestr = "Stokes I signal-to-noise"
                 if levels == "Default":
                     levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(SNRi[stk_cov[0, 0] > 0.0])
-                cont_stk = self.ax_overplot.contour(
-                    SNRi,
-                    levels=levels,
-                    colors="grey",
-                    alpha=0.75,
-                    transform=self.ax_overplot.get_transform(self.wcs_UV),
-                )
+                cont_stk = self.ax_overplot.contour(SNRi, levels=levels, colors="grey", alpha=0.75, transform=self.ax_overplot.get_transform(self.wcs_UV))
             else:  # default to intensity contours
                 disptypestr = "Stokes I"
                 if levels == "Default":
                     levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(stkI[stkI > 0.0]) * self.map_convert
                 cont_stk = self.ax_overplot.contour(
-                    stkI * self.map_convert,
-                    levels=levels,
-                    colors="grey",
-                    alpha=0.75,
-                    transform=self.ax_overplot.get_transform(self.wcs_UV),
+                    stkI * self.map_convert, levels=levels, colors="grey", alpha=0.75, transform=self.ax_overplot.get_transform(self.wcs_UV)
                 )
         # self.ax_overplot.clabel(cont_stk, inline=False, colors="k", fontsize=7)
 
@@ -1861,11 +1579,7 @@ class overplot_pol(align_maps):
             )
             self.ax_overplot.add_artist(pol_sc)
 
-        (self.cr_map,) = self.ax_overplot.plot(
-            *(self.map_wcs.celestial.wcs.crpix - (1.0, 1.0)),
-            "r+",
-            transform=self.ax_overplot.get_transform(self.wcs_UV),
-        )
+        (self.cr_map,) = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix - (1.0, 1.0)), "r+", transform=self.ax_overplot.get_transform(self.wcs_UV))
         (self.cr_other,) = self.ax_overplot.plot(*(self.other_wcs.celestial.wcs.crpix - (1.0, 1.0)), "g+")
 
         if "PHOTPLAM" in list(self.other_header.keys()):
@@ -1885,12 +1599,7 @@ class overplot_pol(align_maps):
             handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=cont_stk.get_edgecolor()[0]))
             disptypestr += " contours"
         self.legend = self.ax_overplot.legend(
-            handles=handles,
-            labels=labels,
-            bbox_to_anchor=(0.0, 1.02, 1.0, 0.102),
-            loc="lower left",
-            mode="expand",
-            borderaxespad=0.0,
+            handles=handles, labels=labels, bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), loc="lower left", mode="expand", borderaxespad=0.0
         )
 
         self.fig_overplot.suptitle(
@@ -1905,29 +1614,10 @@ class overplot_pol(align_maps):
 
         self.fig_overplot.canvas.draw()
 
-    def plot(
-        self,
-        levels=None,
-        P_cut=0.99,
-        SNRi_cut=1.0,
-        step_vec=1,
-        scale_vec=2.0,
-        disptype="i",
-        savename=None,
-        **kwargs,
-    ) -> None:
+    def plot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, step_vec=1, scale_vec=2.0, disptype="i", savename=None, **kwargs) -> None:
         while not self.aligned:
             self.align()
-        self.overplot(
-            levels=levels,
-            P_cut=P_cut,
-            SNRi_cut=SNRi_cut,
-            step_vec=step_vec,
-            scale_vec=scale_vec,
-            disptype=disptype,
-            savename=savename,
-            **kwargs,
-        )
+        self.overplot(levels=levels, P_cut=P_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, scale_vec=scale_vec, disptype=disptype, savename=savename, **kwargs)
         plt.show(block=True)
 
     def add_vector(self, position="center", pol_deg=1.0, pol_ang=0.0, **kwargs):
@@ -1936,10 +1626,7 @@ class overplot_pol(align_maps):
         if isinstance(position, SkyCoord):
             position = self.other_wcs.world_to_pixel(position)
 
-        u, v = (
-            pol_deg * np.cos(np.radians(pol_ang) + np.pi / 2.0),
-            pol_deg * np.sin(np.radians(pol_ang) + np.pi / 2.0),
-        )
+        u, v = (pol_deg * np.cos(np.radians(pol_ang) + np.pi / 2.0), pol_deg * np.sin(np.radians(pol_ang) + np.pi / 2.0))
         for key, value in [
             ["scale", [["scale", 1.0 / self.scale_vec]]],
             ["width", [["width", 0.5 * self.px_scale]]],
@@ -1952,34 +1639,15 @@ class overplot_pol(align_maps):
             except KeyError:
                 for key_i, val_i in value:
                     kwargs[key_i] = val_i
-        handles, labels = (
-            self.legend.legend_handles,
-            [l.get_text() for l in self.legend.texts],
-        )
+        handles, labels = (self.legend.legend_handles, [l.get_text() for l in self.legend.texts])
         new_vec = self.ax_overplot.quiver(
-            *position,
-            u,
-            v,
-            units="xy",
-            angles="uv",
-            scale_units="xy",
-            pivot="mid",
-            headwidth=0.0,
-            headlength=0.0,
-            headaxislength=0.0,
-            **kwargs,
+            *position, u, v, units="xy", angles="uv", scale_units="xy", pivot="mid", headwidth=0.0, headlength=0.0, headaxislength=0.0, **kwargs
         )
         handles.append(FancyArrowPatch((0, 0), (0, 1), arrowstyle="-", fc=new_vec.get_fc(), ec=new_vec.get_ec()))
         labels.append(new_vec.get_label())
         self.legend.remove()
         self.legend = self.ax_overplot.legend(
-            title=self.legend_title,
-            handles=handles,
-            labels=labels,
-            bbox_to_anchor=(0.0, 1.02, 1.0, 0.102),
-            loc="lower left",
-            mode="expand",
-            borderaxespad=0.0,
+            title=self.legend_title, handles=handles, labels=labels, bbox_to_anchor=(0.0, 1.02, 1.0, 0.102), loc="lower left", mode="expand", borderaxespad=0.0
         )
         self.fig_overplot.canvas.draw()
         return new_vec
@@ -1998,17 +1666,7 @@ class align_pol(object):
 
         self.kwargs = kwargs
 
-    def single_plot(
-        self,
-        curr_map,
-        wcs,
-        v_lim=None,
-        ax_lim=None,
-        P_cut=3.0,
-        SNRi_cut=3.0,
-        savename=None,
-        **kwargs,
-    ):
+    def single_plot(self, curr_map, wcs, v_lim=None, ax_lim=None, P_cut=3.0, SNRi_cut=3.0, savename=None, **kwargs):
         # Get data
         stkI = curr_map["I_STOKES"].data
         stk_cov = curr_map["IQU_COV_MATRIX"].data
@@ -2057,10 +1715,7 @@ class align_pol(object):
         else:
             vmin, vmax = v_lim * convert_flux
 
-        for key, value in [
-            ["cmap", [["cmap", "inferno"]]],
-            ["norm", [["vmin", vmin], ["vmax", vmax]]],
-        ]:
+        for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["vmin", vmin], ["vmax", vmax]]]]:
             try:
                 test = kwargs[key]
                 if isinstance(test, LogNorm):
@@ -2070,28 +1725,10 @@ class align_pol(object):
                     kwargs[key_i] = val_i
 
         im = ax.imshow(stkI * convert_flux, aspect="equal", **kwargs)
-        fig.colorbar(
-            im,
-            ax=ax,
-            aspect=50,
-            shrink=0.75,
-            pad=0.025,
-            label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]",
-        )
+        fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
 
         px_size = wcs.wcs.get_cdelt()[0] * 3600.0
-        px_sc = AnchoredSizeBar(
-            ax.transData,
-            1.0 / px_size,
-            "1 arcsec",
-            3,
-            pad=0.5,
-            sep=5,
-            borderpad=0.5,
-            frameon=False,
-            size_vertical=0.005,
-            color="w",
-        )
+        px_sc = AnchoredSizeBar(ax.transData, 1.0 / px_size, "1 arcsec", 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color="w")
         ax.add_artist(px_sc)
 
         north_dir = AnchoredDirectionArrows(
@@ -2116,10 +1753,7 @@ class align_pol(object):
 
         step_vec = 1
         X, Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
-        U, V = (
-            pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0),
-            pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0),
-        )
+        U, V = (pol * np.cos(np.pi / 2.0 + pang * np.pi / 180.0), pol * np.sin(np.pi / 2.0 + pang * np.pi / 180.0))
         ax.quiver(
             X[::step_vec, ::step_vec],
             Y[::step_vec, ::step_vec],
@@ -2137,18 +1771,7 @@ class align_pol(object):
             linewidth=0.75,
             color="w",
         )
-        pol_sc = AnchoredSizeBar(
-            ax.transData,
-            2.0,
-            r"$P$= 100 %",
-            4,
-            pad=0.5,
-            sep=5,
-            borderpad=0.5,
-            frameon=False,
-            size_vertical=0.005,
-            color="w",
-        )
+        pol_sc = AnchoredSizeBar(ax.transData, 2.0, r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color="w")
         ax.add_artist(pol_sc)
 
         if "PHOTPLAM" in list(curr_map[0].header.keys()):
@@ -2183,17 +1806,7 @@ class align_pol(object):
             np.min(
                 [
                     np.min(
-                        curr_map[0].data[
-                            curr_map[0].data
-                            > SNRi_cut
-                            * np.max(
-                                [
-                                    eps * np.ones(curr_map[0].data.shape),
-                                    np.sqrt(curr_map[3].data[0, 0]),
-                                ],
-                                axis=0,
-                            )
-                        ]
+                        curr_map[0].data[curr_map[0].data > SNRi_cut * np.max([eps * np.ones(curr_map[0].data.shape), np.sqrt(curr_map[3].data[0, 0])], axis=0)]
                     )
                     for curr_map in self.other_maps
                 ]
@@ -2202,19 +1815,7 @@ class align_pol(object):
         )
         vmax = np.max(
             [
-                np.max(
-                    curr_map[0].data[
-                        curr_map[0].data
-                        > SNRi_cut
-                        * np.max(
-                            [
-                                eps * np.ones(curr_map[0].data.shape),
-                                np.sqrt(curr_map[3].data[0, 0]),
-                            ],
-                            axis=0,
-                        )
-                    ]
-                )
+                np.max(curr_map[0].data[curr_map[0].data > SNRi_cut * np.max([eps * np.ones(curr_map[0].data.shape), np.sqrt(curr_map[3].data[0, 0])], axis=0)])
                 for curr_map in self.other_maps
             ]
         )
@@ -2224,15 +1825,7 @@ class align_pol(object):
                     vmin,
                     np.min(
                         self.ref_map[0].data[
-                            self.ref_map[0].data
-                            > SNRi_cut
-                            * np.max(
-                                [
-                                    eps * np.ones(self.ref_map[0].data.shape),
-                                    np.sqrt(self.ref_map[3].data[0, 0]),
-                                ],
-                                axis=0,
-                            )
+                            self.ref_map[0].data > SNRi_cut * np.max([eps * np.ones(self.ref_map[0].data.shape), np.sqrt(self.ref_map[3].data[0, 0])], axis=0)
                         ]
                     ),
                 ]
@@ -2244,43 +1837,20 @@ class align_pol(object):
                 vmax,
                 np.max(
                     self.ref_map[0].data[
-                        self.ref_map[0].data
-                        > SNRi_cut
-                        * np.max(
-                            [
-                                eps * np.ones(self.ref_map[0].data.shape),
-                                np.sqrt(self.ref_map[3].data[0, 0]),
-                            ],
-                            axis=0,
-                        )
+                        self.ref_map[0].data > SNRi_cut * np.max([eps * np.ones(self.ref_map[0].data.shape), np.sqrt(self.ref_map[3].data[0, 0])], axis=0)
                     ]
                 ),
             ]
         )
         v_lim = np.array([vmin, vmax])
 
-        fig, ax = self.single_plot(
-            self.ref_map,
-            self.wcs,
-            v_lim=v_lim,
-            P_cut=P_cut,
-            SNRi_cut=SNRi_cut,
-            savename=savename + "_0",
-            **kwargs,
-        )
+        fig, ax = self.single_plot(self.ref_map, self.wcs, v_lim=v_lim, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=savename + "_0", **kwargs)
         x_lim, y_lim = ax.get_xlim(), ax.get_ylim()
         ax_lim = np.array([self.wcs.pixel_to_world(x_lim[i], y_lim[i]) for i in range(len(x_lim))])
 
         for i, curr_map in enumerate(self.other_maps):
             self.single_plot(
-                curr_map,
-                self.wcs_other[i],
-                v_lim=v_lim,
-                ax_lim=ax_lim,
-                P_cut=P_cut,
-                SNRi_cut=SNRi_cut,
-                savename=savename + "_" + str(i + 1),
-                **kwargs,
+                curr_map, self.wcs_other[i], v_lim=v_lim, ax_lim=ax_lim, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=savename + "_" + str(i + 1), **kwargs
             )
 
 
@@ -2347,10 +1917,7 @@ class crop_map(object):
         else:
             kwargs = {**self.kwargs, **kwargs}
 
-        vmin, vmax = (
-            np.min(data[data > 0.0] * convert_flux),
-            np.max(data[data > 0.0] * convert_flux),
-        )
+        vmin, vmax = (np.min(data[data > 0.0] * convert_flux), np.max(data[data > 0.0] * convert_flux))
         for key, value in [
             ["cmap", [["cmap", "inferno"]]],
             ["origin", [["origin", "lower"]]],
@@ -2563,15 +2130,9 @@ class crop_Stokes(crop_map):
             if dataset.header["FILENAME"][-4:] != "crop":
                 dataset.header["FILENAME"] += "_crop"
             dataset.header["P_int"] = (P_diluted, "Integrated polarization degree")
-            dataset.header["sP_int"] = (
-                np.ceil(P_diluted_err * 1000.0) / 1000.0,
-                "Integrated polarization degree error",
-            )
+            dataset.header["sP_int"] = (np.ceil(P_diluted_err * 1000.0) / 1000.0, "Integrated polarization degree error")
             dataset.header["PA_int"] = (PA_diluted, "Integrated polarization angle")
-            dataset.header["sPA_int"] = (
-                np.ceil(PA_diluted_err * 10.0) / 10.0,
-                "Integrated polarization angle error",
-            )
+            dataset.header["sPA_int"] = (np.ceil(PA_diluted_err * 10.0) / 10.0, "Integrated polarization angle error")
         self.fig.canvas.draw_idle()
 
     @property
@@ -2600,14 +2161,7 @@ class image_lasso_selector(object):
             self.ax = ax
             self.mask_alpha = 0.1
             self.embedded = True
-        self.displayed = self.ax.imshow(
-            self.img,
-            vmin=self.vmin,
-            vmax=self.vmax,
-            aspect="equal",
-            cmap="inferno",
-            alpha=self.mask_alpha,
-        )
+        self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect="equal", cmap="inferno", alpha=self.mask_alpha)
         plt.ion()
 
         lineprops = {"color": "grey", "linewidth": 1, "alpha": 0.8}
@@ -2636,14 +2190,7 @@ class image_lasso_selector(object):
 
     def update_mask(self):
         self.displayed.remove()
-        self.displayed = self.ax.imshow(
-            self.img,
-            vmin=self.vmin,
-            vmax=self.vmax,
-            aspect="equal",
-            cmap="inferno",
-            alpha=self.mask_alpha,
-        )
+        self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect="equal", cmap="inferno", alpha=self.mask_alpha)
         array = self.displayed.get_array().data
 
         self.mask = np.zeros(self.img.shape[:2], dtype=bool)
@@ -2659,16 +2206,7 @@ class image_lasso_selector(object):
 
 
 class slit(object):
-    def __init__(
-        self,
-        img,
-        cdelt=np.array([1.0, 1.0]),
-        width=1.0,
-        height=2.0,
-        angle=0.0,
-        fig=None,
-        ax=None,
-    ):
+    def __init__(self, img, cdelt=np.array([1.0, 1.0]), width=1.0, height=2.0, angle=0.0, fig=None, ax=None):
         """
         img must have shape (X, Y)
         """
@@ -2689,14 +2227,7 @@ class slit(object):
             self.mask_alpha = 0.1
             self.embedded = True
 
-        self.displayed = self.ax.imshow(
-            self.img,
-            vmin=self.vmin,
-            vmax=self.vmax,
-            aspect="equal",
-            cmap="inferno",
-            alpha=self.mask_alpha,
-        )
+        self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect="equal", cmap="inferno", alpha=self.mask_alpha)
         plt.ion()
 
         xx, yy = np.indices(self.img.shape)
@@ -2710,15 +2241,7 @@ class slit(object):
         self.angle = angle
 
         self.rect_center = (self.x0, self.y0) - np.dot(rot2D(self.angle), (self.width / 2, self.height / 2))
-        self.rect = Rectangle(
-            self.rect_center,
-            self.width,
-            self.height,
-            angle=self.angle,
-            alpha=0.8,
-            ec="grey",
-            fc="none",
-        )
+        self.rect = Rectangle(self.rect_center, self.width, self.height, angle=self.angle, alpha=0.8, ec="grey", fc="none")
         self.ax.add_patch(self.rect)
 
         self.fig.canvas.mpl_connect("button_press_event", self.on_press)
@@ -2778,14 +2301,7 @@ class slit(object):
                 self.displayed.remove()
             except ValueError:
                 return
-        self.displayed = self.ax.imshow(
-            self.img,
-            vmin=self.vmin,
-            vmax=self.vmax,
-            aspect="equal",
-            cmap="inferno",
-            alpha=self.mask_alpha,
-        )
+        self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect="equal", cmap="inferno", alpha=self.mask_alpha)
         array = self.displayed.get_array().data
 
         self.mask = np.zeros(array.shape, dtype=bool)
@@ -2823,14 +2339,7 @@ class aperture(object):
             self.mask_alpha = 0.1
             self.embedded = True
 
-        self.displayed = self.ax.imshow(
-            self.img,
-            vmin=self.vmin,
-            vmax=self.vmax,
-            aspect="equal",
-            cmap="inferno",
-            alpha=self.mask_alpha,
-        )
+        self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect="equal", cmap="inferno", alpha=self.mask_alpha)
         plt.ion()
 
         xx, yy = np.indices(self.img.shape)
@@ -2891,14 +2400,7 @@ class aperture(object):
                 self.displayed.remove()
             except ValueError:
                 return
-        self.displayed = self.ax.imshow(
-            self.img,
-            vmin=self.vmin,
-            vmax=self.vmax,
-            aspect="equal",
-            cmap="inferno",
-            alpha=self.mask_alpha,
-        )
+        self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect="equal", cmap="inferno", alpha=self.mask_alpha)
         array = self.displayed.get_array().data
 
         yy, xx = np.indices(self.img.shape[:2])
@@ -2919,17 +2421,7 @@ class pol_map(object):
     Class to interactively study polarization maps.
     """
 
-    def __init__(
-        self,
-        Stokes,
-        P_cut=0.99,
-        SNRi_cut=1.0,
-        step_vec=1,
-        scale_vec=3.0,
-        flux_lim=None,
-        selection=None,
-        pa_err=False,
-    ):
+    def __init__(self, Stokes, P_cut=0.99, SNRi_cut=1.0, step_vec=1, scale_vec=3.0, flux_lim=None, selection=None, pa_err=False):
         if isinstance(Stokes, str):
             Stokes = fits.open(Stokes)
         self.Stokes = deepcopy(Stokes)
@@ -2972,22 +2464,8 @@ class pol_map(object):
         ax_snr_conf = self.fig.add_axes([0.115, 0.020, 0.05, 0.02])
         SNRi_max = np.max(self.I[self.IQU_cov[0, 0] > 0.0] / np.sqrt(self.IQU_cov[0, 0][self.IQU_cov[0, 0] > 0.0]))
         SNRp_max = np.max(self.P[self.s_P > 0.0] / self.s_P[self.s_P > 0.0])
-        s_I_cut = Slider(
-            ax_I_cut,
-            r"$SNR^{I}_{cut}$",
-            1.0,
-            int(SNRi_max * 0.95),
-            valstep=1,
-            valinit=self.SNRi_cut,
-        )
-        self.s_P_cut = Slider(
-            self.ax_P_cut,
-            r"$Conf^{P}_{cut}$",
-            0.50,
-            1.0,
-            valstep=[0.500, 0.900, 0.990, 0.999],
-            valinit=self.P_cut if P_cut <= 1.0 else 0.99,
-        )
+        s_I_cut = Slider(ax_I_cut, r"$SNR^{I}_{cut}$", 1.0, int(SNRi_max * 0.95), valstep=1, valinit=self.SNRi_cut)
+        self.s_P_cut = Slider(self.ax_P_cut, r"$Conf^{P}_{cut}$", 0.50, 1.0, valstep=[0.500, 0.900, 0.990, 0.999], valinit=self.P_cut if P_cut <= 1.0 else 0.99)
         s_vec_sc = Slider(ax_vec_sc, r"Vec scale", 0.0, 10.0, valstep=1, valinit=self.scale_vec)
         b_snr_reset = Button(ax_snr_reset, "Reset")
         b_snr_reset.label.set_fontsize(8)
@@ -3031,12 +2509,7 @@ class pol_map(object):
                 self.snr_conf = 1
                 b_snr_conf.label.set_text("SNR")
                 self.s_P_cut = Slider(
-                    self.ax_P_cut,
-                    r"$Conf^{P}_{cut}$",
-                    0.50,
-                    1.0,
-                    valstep=[0.500, 0.900, 0.990, 0.999],
-                    valinit=self.P_cut if P_cut <= 1.0 else 0.99,
+                    self.ax_P_cut, r"$Conf^{P}_{cut}$", 0.50, 1.0, valstep=[0.500, 0.900, 0.990, 0.999], valinit=self.P_cut if P_cut <= 1.0 else 0.99
                 )
             self.s_P_cut.on_changed(update_snrp)
             update_snrp(self.s_P_cut.val)
@@ -3100,14 +2573,7 @@ class pol_map(object):
         b_aper.label.set_fontsize(8)
         b_aper_reset = Button(ax_aper_reset, "Reset")
         b_aper_reset.label.set_fontsize(8)
-        s_aper_radius = Slider(
-            ax_aper_radius,
-            r"$R_{aper}$",
-            np.ceil(self.wcs.wcs.cdelt.max() / 1.33 * 3.6e5) / 1e2,
-            3.5,
-            valstep=1e-2,
-            valinit=1.0,
-        )
+        s_aper_radius = Slider(ax_aper_radius, r"$R_{aper}$", np.ceil(self.wcs.wcs.cdelt.max() / 1.33 * 3.6e5) / 1e2, 3.5, valstep=1e-2, valinit=1.0)
 
         def select_aperture(event):
             if self.data is None:
@@ -3124,13 +2590,7 @@ class pol_map(object):
             else:
                 self.selected = True
                 self.region = None
-                self.select_instance = aperture(
-                    self.data,
-                    fig=self.fig,
-                    ax=self.ax,
-                    cdelt=self.wcs.wcs.cdelt,
-                    radius=s_aper_radius.val,
-                )
+                self.select_instance = aperture(self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, radius=s_aper_radius.val)
                 self.select_instance.circ.set_visible(True)
 
             self.fig.canvas.draw_idle()
@@ -3141,22 +2601,10 @@ class pol_map(object):
                     self.select_instance.update_radius(val)
                 else:
                     self.selected = True
-                    self.select_instance = aperture(
-                        self.data,
-                        fig=self.fig,
-                        ax=self.ax,
-                        cdelt=self.wcs.wcs.cdelt,
-                        radius=val,
-                    )
+                    self.select_instance = aperture(self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, radius=val)
             else:
                 self.selected = True
-                self.select_instance = aperture(
-                    self.data,
-                    fig=self.fig,
-                    ax=self.ax,
-                    cdelt=self.wcs.wcs.cdelt,
-                    radius=val,
-                )
+                self.select_instance = aperture(self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, radius=val)
             self.fig.canvas.draw_idle()
 
         def reset_aperture(event):
@@ -3180,22 +2628,8 @@ class pol_map(object):
         b_slit.label.set_fontsize(8)
         b_slit_reset = Button(ax_slit_reset, "Reset")
         b_slit_reset.label.set_fontsize(8)
-        s_slit_width = Slider(
-            ax_slit_width,
-            r"$W_{slit}$",
-            np.ceil(self.wcs.wcs.cdelt.max() / 1.33 * 3.6e5) / 1e2,
-            7.0,
-            valstep=1e-2,
-            valinit=1.0,
-        )
-        s_slit_height = Slider(
-            ax_slit_height,
-            r"$H_{slit}$",
-            np.ceil(self.wcs.wcs.cdelt.max() / 1.33 * 3.6e5) / 1e2,
-            7.0,
-            valstep=1e-2,
-            valinit=1.0,
-        )
+        s_slit_width = Slider(ax_slit_width, r"$W_{slit}$", np.ceil(self.wcs.wcs.cdelt.max() / 1.33 * 3.6e5) / 1e2, 7.0, valstep=1e-2, valinit=1.0)
+        s_slit_height = Slider(ax_slit_height, r"$H_{slit}$", np.ceil(self.wcs.wcs.cdelt.max() / 1.33 * 3.6e5) / 1e2, 7.0, valstep=1e-2, valinit=1.0)
         s_slit_angle = Slider(ax_slit_angle, r"$\theta_{slit}$", 0.0, 90.0, valstep=1.0, valinit=0.0)
 
         def select_slit(event):
@@ -3214,13 +2648,7 @@ class pol_map(object):
                 self.selected = True
                 self.region = None
                 self.select_instance = slit(
-                    self.data,
-                    fig=self.fig,
-                    ax=self.ax,
-                    cdelt=self.wcs.wcs.cdelt,
-                    width=s_slit_width.val,
-                    height=s_slit_height.val,
-                    angle=s_slit_angle.val,
+                    self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=s_slit_width.val, height=s_slit_height.val, angle=s_slit_angle.val
                 )
                 self.select_instance.rect.set_visible(True)
 
@@ -3233,24 +2661,12 @@ class pol_map(object):
                 else:
                     self.selected = True
                     self.select_instance = slit(
-                        self.data,
-                        fig=self.fig,
-                        ax=self.ax,
-                        cdelt=self.wcs.wcs.cdelt,
-                        width=val,
-                        height=s_slit_height.val,
-                        angle=s_slit_angle.val,
+                        self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=val, height=s_slit_height.val, angle=s_slit_angle.val
                     )
             else:
                 self.selected = True
                 self.select_instance = slit(
-                    self.data,
-                    fig=self.fig,
-                    ax=self.ax,
-                    cdelt=self.wcs.wcs.cdelt,
-                    width=val,
-                    height=s_slit_height.val,
-                    angle=s_slit_angle.val,
+                    self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=val, height=s_slit_height.val, angle=s_slit_angle.val
                 )
             self.fig.canvas.draw_idle()
 
@@ -3261,24 +2677,12 @@ class pol_map(object):
                 else:
                     self.selected = True
                     self.select_instance = slit(
-                        self.data,
-                        fig=self.fig,
-                        ax=self.ax,
-                        cdelt=self.wcs.wcs.cdelt,
-                        width=s_slit_width.val,
-                        height=val,
-                        angle=s_slit_angle.val,
+                        self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=s_slit_width.val, height=val, angle=s_slit_angle.val
                     )
             else:
                 self.selected = True
                 self.select_instance = slit(
-                    self.data,
-                    fig=self.fig,
-                    ax=self.ax,
-                    cdelt=self.wcs.wcs.cdelt,
-                    width=s_slit_width.val,
-                    height=val,
-                    angle=s_slit_angle.val,
+                    self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=s_slit_width.val, height=val, angle=s_slit_angle.val
                 )
             self.fig.canvas.draw_idle()
 
@@ -3289,24 +2693,12 @@ class pol_map(object):
                 else:
                     self.selected = True
                     self.select_instance = slit(
-                        self.data,
-                        fig=self.fig,
-                        ax=self.ax,
-                        cdelt=self.wcs.wcs.cdelt,
-                        width=s_slit_width.val,
-                        height=s_slit_height.val,
-                        angle=val,
+                        self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=s_slit_width.val, height=s_slit_height.val, angle=val
                     )
             else:
                 self.selected = True
                 self.select_instance = slit(
-                    self.data,
-                    fig=self.fig,
-                    ax=self.ax,
-                    cdelt=self.wcs.wcs.cdelt,
-                    width=s_slit_width.val,
-                    height=s_slit_height.val,
-                    angle=val,
+                    self.data, fig=self.fig, ax=self.ax, cdelt=self.wcs.wcs.cdelt, width=s_slit_width.val, height=s_slit_height.val, angle=val
                 )
             self.fig.canvas.draw_idle()
 
@@ -3391,11 +2783,7 @@ class pol_map(object):
         def submit_save(expression):
             ax_text_save.set(visible=False)
             if expression != "":
-                save_fig, save_ax = plt.subplots(
-                    figsize=(12, 10),
-                    layout="constrained",
-                    subplot_kw=dict(projection=self.wcs),
-                )
+                save_fig, save_ax = plt.subplots(figsize=(12, 10), layout="constrained", subplot_kw=dict(projection=self.wcs))
                 self.ax_cosmetics(ax=save_ax)
                 self.display(fig=save_fig, ax=save_ax)
                 self.pol_vector(fig=save_fig, ax=save_ax)
@@ -3433,10 +2821,7 @@ class pol_map(object):
             center = (shape / 2).astype(int)
             cdelt_arcsec = self.wcs.wcs.cdelt * 3600
             xx, yy = np.indices(shape)
-            x, y = (
-                (xx - center[0]) * cdelt_arcsec[0],
-                (yy - center[1]) * cdelt_arcsec[1],
-            )
+            x, y = ((xx - center[0]) * cdelt_arcsec[0], (yy - center[1]) * cdelt_arcsec[1])
 
             P, PA = np.zeros(shape), np.zeros(shape)
             P[self.cut] = self.P[self.cut]
@@ -3445,15 +2830,7 @@ class pol_map(object):
             for i in range(shape[0]):
                 for j in range(shape[1]):
                     dump_list.append(
-                        [
-                            x[i, j],
-                            y[i, j],
-                            self.I[i, j] * self.map_convert,
-                            self.Q[i, j] * self.map_convert,
-                            self.U[i, j] * self.map_convert,
-                            P[i, j],
-                            PA[i, j],
-                        ]
+                        [x[i, j], y[i, j], self.I[i, j] * self.map_convert, self.Q[i, j] * self.map_convert, self.U[i, j] * self.map_convert, P[i, j], PA[i, j]]
                     )
             self.data_dump = np.array(dump_list)
 
@@ -3606,10 +2983,7 @@ class pol_map(object):
     @property
     def cut(self):
         s_I = np.sqrt(self.IQU_cov[0, 0])
-        SNRp_mask, SNRi_mask = (
-            np.zeros(self.P.shape).astype(bool),
-            np.zeros(self.I.shape).astype(bool),
-        )
+        SNRp_mask, SNRi_mask = (np.zeros(self.P.shape).astype(bool), np.zeros(self.I.shape).astype(bool))
         SNRi_mask[s_I > 0.0] = self.I[s_I > 0.0] / s_I[s_I > 0.0] > self.SNRi
         if self.SNRp >= 1.0:
             SNRp_mask[self.s_P > 0.0] = self.P[self.s_P > 0.0] / self.s_P[self.s_P > 0.0] > self.SNRp
@@ -3654,17 +3028,7 @@ class pol_map(object):
         else:
             scale_vec = 2.0
         self.pol_sc = AnchoredSizeBar(
-            ax.transData,
-            scale_vec,
-            r"$P$= 100%",
-            4,
-            pad=0.5,
-            sep=5,
-            borderpad=0.5,
-            frameon=False,
-            size_vertical=0.005,
-            color="white",
-            fontproperties=fontprops,
+            ax.transData, scale_vec, r"$P$= 100%", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color="white", fontproperties=fontprops
         )
         ax.add_artist(self.pol_sc)
         if hasattr(self, "north_dir"):
@@ -3698,10 +3062,7 @@ class pol_map(object):
         if self.display_selection.lower() in ["total_flux"]:
             self.data = self.I * self.map_convert
             if flux_lim is None:
-                vmin, vmax = (
-                    1.0 / 2.0 * np.median(self.data[self.data > 0.0]),
-                    np.max(self.data[self.data > 0.0]),
-                )
+                vmin, vmax = (1.0 / 2.0 * np.median(self.data[self.data > 0.0]), np.max(self.data[self.data > 0.0]))
             else:
                 vmin, vmax = flux_lim
             norm = LogNorm(vmin, vmax)
@@ -3709,10 +3070,7 @@ class pol_map(object):
         elif self.display_selection.lower() in ["pol_flux"]:
             self.data = self.I * self.map_convert * self.P
             if flux_lim is None:
-                vmin, vmax = (
-                    1.0 / 2.0 * np.median(self.I[self.I > 0.0] * self.map_convert),
-                    np.max(self.I[self.I > 0.0] * self.map_convert),
-                )
+                vmin, vmax = (1.0 / 2.0 * np.median(self.I[self.I > 0.0] * self.map_convert), np.max(self.I[self.I > 0.0] * self.map_convert))
             else:
                 vmin, vmax = flux_lim
             norm = LogNorm(vmin, vmax)
@@ -3778,10 +3136,7 @@ class pol_map(object):
             scale_vec = 2.0
             P_cut[self.cut] = 1.0 / scale_vec
         X, Y = np.meshgrid(np.arange(self.I.shape[1]), np.arange(self.I.shape[0]))
-        XY_U, XY_V = (
-            P_cut * np.cos(np.pi / 2.0 + self.PA * np.pi / 180.0),
-            P_cut * np.sin(np.pi / 2.0 + self.PA * np.pi / 180.0),
-        )
+        XY_U, XY_V = (P_cut * np.cos(np.pi / 2.0 + self.PA * np.pi / 180.0), P_cut * np.sin(np.pi / 2.0 + self.PA * np.pi / 180.0))
 
         if fig is None:
             fig = self.fig
@@ -3985,12 +3340,7 @@ class pol_map(object):
             IU_reg_err = np.sqrt(np.sum(s_IU[self.region] ** 2))
             QU_reg_err = np.sqrt(np.sum(s_QU[self.region] ** 2))
 
-            conf = PCconf(
-                QN=Q_reg / I_reg,
-                QN_ERR=Q_reg_err / I_reg,
-                UN=U_reg / I_reg,
-                UN_ERR=U_reg_err / I_reg,
-            )
+            conf = PCconf(QN=Q_reg / I_reg, QN_ERR=Q_reg_err / I_reg, UN=U_reg / I_reg, UN_ERR=U_reg_err / I_reg)
             if 1.0 - conf > 1e-3:
                 str_conf = "\n" + r"Confidence: {0:.2f} %".format(conf * 100.0)
 
@@ -4050,8 +3400,7 @@ class pol_map(object):
                 self.an_int.remove()
             self.str_int = (
                 r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(
-                    self.pivot_wav,
-                    sci_not(I_reg * self.map_convert, I_reg_err * self.map_convert, 2),
+                    self.pivot_wav, sci_not(I_reg * self.map_convert, I_reg_err * self.map_convert, 2)
                 )
                 + "\n"
                 + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg * 100.0, np.ceil(P_reg_err * 1000.0) / 10.0)
@@ -4081,19 +3430,13 @@ class pol_map(object):
                 horizontalalignment="left",
             )
             if self.region is not None:
-                self.cont = ax.contour(
-                    self.region.astype(float),
-                    levels=[0.5],
-                    colors="white",
-                    linewidths=0.8,
-                )
+                self.cont = ax.contour(self.region.astype(float), levels=[0.5], colors="white", linewidths=0.8)
             fig.canvas.draw_idle()
             return self.an_int
         else:
             str_int = (
                 r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(
-                    self.pivot_wav,
-                    sci_not(I_reg * self.map_convert, I_reg_err * self.map_convert, 2),
+                    self.pivot_wav, sci_not(I_reg * self.map_convert, I_reg_err * self.map_convert, 2)
                 )
                 + "\n"
                 + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg * 100.0, np.ceil(P_reg_err * 1000.0) / 10.0)
@@ -4123,12 +3466,7 @@ class pol_map(object):
                 horizontalalignment="left",
             )
             if self.region is not None:
-                ax.contour(
-                    self.region.astype(float),
-                    levels=[0.5],
-                    colors="white",
-                    linewidths=0.8,
-                )
+                ax.contour(self.region.astype(float), levels=[0.5], colors="white", linewidths=0.8)
             fig.canvas.draw_idle()
 
 
@@ -4136,84 +3474,21 @@ if __name__ == "__main__":
     import argparse
 
     parser = argparse.ArgumentParser(description="Interactively plot the pipeline products")
+    parser.add_argument("-f", "--file", metavar="path", required=False, help="The full or relative path to the data product", type=str, default=None)
+    parser.add_argument("-p", "--pcut", metavar="p_cut", required=False, help="The cut in signal-to-noise for the polarization degree", type=float, default=3.0)
+    parser.add_argument("-i", "--snri", metavar="snri_cut", required=False, help="The cut in signal-to-noise for the intensity", type=float, default=3.0)
     parser.add_argument(
-        "-f",
-        "--file",
-        metavar="path",
-        required=False,
-        help="The full or relative path to the data product",
-        type=str,
-        default=None,
+        "-st", "--step-vec", metavar="step_vec", required=False, help="Quantity of vectors to be shown, 1 is all, 2 is every other, etc.", type=int, default=1
     )
     parser.add_argument(
-        "-p",
-        "--pcut",
-        metavar="p_cut",
-        required=False,
-        help="The cut in signal-to-noise for the polarization degree",
-        type=float,
-        default=3.0,
+        "-sc", "--scale-vec", metavar="scale_vec", required=False, help="Size of the 100%% polarization vector in pixel units", type=float, default=3.0
     )
+    parser.add_argument("-pa", "--pang-err", action="store_true", required=False, help="Whether the polarization angle uncertainties should be displayed")
+    parser.add_argument("-l", "--lim", metavar="flux_lim", nargs=2, required=False, help="Limits for the intensity map", type=float, default=None)
     parser.add_argument(
-        "-i",
-        "--snri",
-        metavar="snri_cut",
-        required=False,
-        help="The cut in signal-to-noise for the intensity",
-        type=float,
-        default=3.0,
-    )
-    parser.add_argument(
-        "-st",
-        "--step-vec",
-        metavar="step_vec",
-        required=False,
-        help="Quantity of vectors to be shown, 1 is all, 2 is every other, etc.",
-        type=int,
-        default=1,
-    )
-    parser.add_argument(
-        "-sc",
-        "--scale-vec",
-        metavar="scale_vec",
-        required=False,
-        help="Size of the 100%% polarization vector in pixel units",
-        type=float,
-        default=3.0,
-    )
-    parser.add_argument(
-        "-pa",
-        "--pang-err",
-        action="store_true",
-        required=False,
-        help="Whether the polarization angle uncertainties should be displayed",
-    )
-    parser.add_argument(
-        "-l",
-        "--lim",
-        metavar="flux_lim",
-        nargs=2,
-        required=False,
-        help="Limits for the intensity map",
-        type=float,
-        default=None,
-    )
-    parser.add_argument(
-        "-pdf",
-        "--static-pdf",
-        metavar="static_pdf",
-        required=False,
-        help="Whether the analysis tool or the static pdfs should be outputed",
-        default=None,
-    )
-    parser.add_argument(
-        "-t",
-        "--type",
-        metavar="type",
-        required=False,
-        help="Type of plot to be be outputed",
-        default=None,
+        "-pdf", "--static-pdf", metavar="static_pdf", required=False, help="Whether the analysis tool or the static pdfs should be outputed", default=None
     )
+    parser.add_argument("-t", "--type", metavar="type", required=False, help="Type of plot to be be outputed", default=None)
     args = parser.parse_args()
 
     if args.file is not None:
diff --git a/package/lib/query.py b/package/lib/query.py
index edfbf2f..3a0aca4 100755
--- a/package/lib/query.py
+++ b/package/lib/query.py
@@ -173,11 +173,7 @@ def get_product_list(target=None, proposal_id=None, instrument="foc"):
 
     observations = Observations.query_criteria(obs_id=list(results["Dataset"][b]))
     products = Observations.filter_products(
-        Observations.get_product_list(observations),
-        productType=["SCIENCE"],
-        dataproduct_type=dataproduct_type,
-        calib_level=[2],
-        description=description,
+        Observations.get_product_list(observations), productType=["SCIENCE"], dataproduct_type=dataproduct_type, calib_level=[2], description=description
     )
 
     products["proposal_id"] = Column(products["proposal_id"], dtype="U35")
diff --git a/package/lib/reduction.py b/package/lib/reduction.py
index b663c3d..89d52b5 100755
--- a/package/lib/reduction.py
+++ b/package/lib/reduction.py
@@ -191,8 +191,8 @@ def bin_ndarray(ndarray, new_shape, operation="sum"):
 
     Example
     -------
-    >>> m = np.arange(0,100,1).reshape((10,10))
-    >>> n = bin_ndarray(m, new_shape=(5,5), operation='sum')
+    >>> m = np.arange(0, 100, 1).reshape((10, 10))
+    >>> n = bin_ndarray(m, new_shape=(5, 5), operation="sum")
     >>> print(n)
 
     [[ 22  30  38  46  54]
@@ -279,9 +279,7 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu
     if null_val is None:
         null_val = [1.00 * error.mean() for error in error_array]
     elif type(null_val) is float:
-        null_val = [
-            null_val,
-        ] * error_array.shape[0]
+        null_val = [null_val] * error_array.shape[0]
 
     vertex = np.zeros((data_array.shape[0], 4), dtype=int)
     for i, image in enumerate(data_array):  # Get vertex of the rectangular convex hull of each image
@@ -348,10 +346,7 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu
                 headers,
                 vmin=convert_flux * data_array[data_array > 0.0].mean() / 5.0,
                 vmax=convert_flux * data_array[data_array > 0.0].max(),
-                rectangle=[
-                    rectangle,
-                ]
-                * len(headers),
+                rectangle=[rectangle] * len(headers),
                 savename=savename + "_crop_region",
                 plots_folder=plots_folder,
             )
@@ -629,12 +624,7 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio
 
         # Compute binning ratio
         if scale.lower() in ["px", "pixel"]:
-            Dxy_arr[i] = np.array(
-                [
-                    pxsize,
-                ]
-                * 2
-            )
+            Dxy_arr[i] = np.array([pxsize] * 2)
             scale = "px"
         elif scale.lower() in ["arcsec", "arcseconds"]:
             Dxy_arr[i] = np.array(pxsize / np.abs(w.wcs.cdelt) / 3600.0)
@@ -938,12 +928,7 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1.5, scale="pi
                 dist_rc = np.where(data_mask, np.sqrt((r - xx) ** 2 + (c - yy) ** 2), fmax)
                 # Catch expected "OverflowWarning" as we overflow values that are not in the image
                 with warnings.catch_warnings(record=True) as w:
-                    g_rc = np.array(
-                        [
-                            np.exp(-0.5 * (dist_rc / stdev) ** 2) / (2.0 * np.pi * stdev**2),
-                        ]
-                        * data_array.shape[0]
-                    )
+                    g_rc = np.array([np.exp(-0.5 * (dist_rc / stdev) ** 2) / (2.0 * np.pi * stdev**2)] * data_array.shape[0])
                     # Apply weighted combination
                     smoothed[r, c] = np.where(data_mask[r, c], np.sum(data_array * weight * g_rc) / np.sum(weight * g_rc), data_array.mean(axis=0)[r, c])
                     error[r, c] = np.where(
@@ -1438,9 +1423,7 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale
         all_Q_stokes = np.zeros((np.unique(rotate).size, data_array.shape[1], data_array.shape[2]))
         all_U_stokes = np.zeros((np.unique(rotate).size, data_array.shape[1], data_array.shape[2]))
         all_Stokes_cov = np.zeros((np.unique(rotate).size, 3, 3, data_array.shape[1], data_array.shape[2]))
-        all_header_stokes = [
-            {},
-        ] * np.unique(rotate).size
+        all_header_stokes = [{}] * np.unique(rotate).size
 
         for i, rot in enumerate(np.unique(rotate)):
             rot_mask = rotate == rot
diff --git a/package/src/emission_center.py b/package/src/emission_center.py
index 097362c..6e3b892 100755
--- a/package/src/emission_center.py
+++ b/package/src/emission_center.py
@@ -23,12 +23,7 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None):
     stkI = Stokes["I_STOKES"].data
     QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan)
     for sflux, nflux in zip(
-        [
-            Stokes["Q_STOKES"].data,
-            Stokes["U_STOKES"].data,
-            np.sqrt(Stokes["IQU_COV_MATRIX"].data[1, 1]),
-            np.sqrt(Stokes["IQU_COV_MATRIX"].data[2, 2]),
-        ],
+        [Stokes["Q_STOKES"].data, Stokes["U_STOKES"].data, np.sqrt(Stokes["IQU_COV_MATRIX"].data[1, 1]), np.sqrt(Stokes["IQU_COV_MATRIX"].data[2, 2])],
         [QN, UN, QN_ERR, UN_ERR],
     ):
         nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0]
@@ -49,25 +44,10 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None):
         target = Stokes[0].header["TARGNAME"]
 
     fig = figure(figsize=(8, 8.5), layout="constrained")
-    fig, ax = polarization_map(
-        Stokes,
-        P_cut=P_cut,
-        step_vec=1,
-        scale_vec=5,
-        display=display,
-        fig=fig,
-        width=0.33,
-        linewidth=0.5,
-    )
+    fig, ax = polarization_map(Stokes, P_cut=P_cut, step_vec=1, scale_vec=5, display=display, fig=fig, width=0.33, linewidth=0.5)
 
     ax.plot(*Stokescenter, marker="+", color="k", lw=3)
-    ax.plot(
-        *Stokescenter,
-        marker="+",
-        color="red",
-        lw=1.5,
-        label="Best confidence for center: {0}".format(Stokespos.to_string("hmsdms")),
-    )
+    ax.plot(*Stokescenter, marker="+", color="red", lw=1.5, label="Best confidence for center: {0}".format(Stokespos.to_string("hmsdms")))
     ax.contour(Stokescentconf, [0.01], colors="k", linewidths=3)
     confcentcont = ax.contour(Stokescentconf, [0.01], colors="red")
     # confcont = ax.contour(Stokesconf, [0.9905], colors="r")
@@ -82,14 +62,7 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None):
     # handles.append(Rectangle((0, 0), 1, 1, fill=False, ls="--", ec=snr3cont.get_edgecolor()[0]))
     # labels.append(r"$SNR_P \geq$ 4  contour")
     # handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=snr4cont.get_edgecolor()[0]))
-    ax.legend(
-        handles=handles,
-        labels=labels,
-        bbox_to_anchor=(-0.05, -0.02, 1.10, 0.01),
-        loc="upper left",
-        mode="expand",
-        borderaxespad=0.0,
-    )
+    ax.legend(handles=handles, labels=labels, bbox_to_anchor=(-0.05, -0.02, 1.10, 0.01), loc="upper left", mode="expand", borderaxespad=0.0)
 
     if output_dir is not None:
         filename = pathjoin(output_dir, "%s_center.pdf" % target)
@@ -103,57 +76,11 @@ if __name__ == "__main__":
     import argparse
 
     parser = argparse.ArgumentParser(description="Look for the center of emission for a given reduced observation")
-    parser.add_argument(
-        "-t",
-        "--target",
-        metavar="targetname",
-        required=False,
-        help="the name of the target",
-        type=str,
-        default=None,
-    )
-    parser.add_argument(
-        "-f",
-        "--file",
-        metavar="path",
-        required=False,
-        help="The full or relative path to the data product",
-        type=str,
-        default=None,
-    )
-    parser.add_argument(
-        "-c",
-        "--pcut",
-        metavar="pcut",
-        required=False,
-        help="The polarization cut for the data mask",
-        type=float,
-        default=0.99,
-    )
-    parser.add_argument(
-        "-d",
-        "--display",
-        metavar="display",
-        required=False,
-        help="The map on which to display info",
-        type=str,
-        default="pf",
-    )
-    parser.add_argument(
-        "-o",
-        "--output_dir",
-        metavar="directory_path",
-        required=False,
-        help="output directory path for the plots",
-        type=str,
-        default="./data",
-    )
+    parser.add_argument("-t", "--target", metavar="targetname", required=False, help="the name of the target", type=str, default=None)
+    parser.add_argument("-f", "--file", metavar="path", required=False, help="The full or relative path to the data product", type=str, default=None)
+    parser.add_argument("-c", "--pcut", metavar="pcut", required=False, help="The polarization cut for the data mask", type=float, default=0.99)
+    parser.add_argument("-d", "--display", metavar="display", required=False, help="The map on which to display info", type=str, default="pf")
+    parser.add_argument("-o", "--output_dir", metavar="directory_path", required=False, help="output directory path for the plots", type=str, default="./data")
     args = parser.parse_args()
-    exitcode = main(
-        infile=args.file,
-        P_cut=args.pcut,
-        target=args.target,
-        display=args.display,
-        output_dir=args.output_dir,
-    )
+    exitcode = main(infile=args.file, P_cut=args.pcut, target=args.target, display=args.display, output_dir=args.output_dir)
     print("Written to: ", exitcode)
diff --git a/package/src/overplot_NGC1068.py b/package/src/overplot_NGC1068.py
index 8642458..3f22923 100755
--- a/package/src/overplot_NGC1068.py
+++ b/package/src/overplot_NGC1068.py
@@ -15,17 +15,7 @@ Radio = fits.open("./data/NGC1068/MERLIN-VLA/Combined_crop.fits")
 
 levels = np.logspace(-0.5, 1.99, 7) / 100.0 * Stokes_UV[0].data.max() * Stokes_UV[0].header["photflam"]
 A = overplot_pol(Stokes_UV, Radio, norm=LogNorm())
-A.plot(
-    levels=levels,
-    P_cut=0.99,
-    SNRi_cut=1.0,
-    scale_vec=3,
-    step_vec=1,
-    norm=LogNorm(5e-5, 1e-1),
-    cmap="inferno_r",
-    width=0.8,
-    linewidth=1.2,
-)
+A.plot(levels=levels, P_cut=0.99, SNRi_cut=1.0, scale_vec=3, step_vec=1, norm=LogNorm(5e-5, 1e-1), cmap="inferno_r", width=0.8, linewidth=1.2)
 A.add_vector(
     A.other_wcs.celestial.wcs.crpix - (1.0, 1.0),
     pol_deg=0.124,
diff --git a/package/src/readfos.py b/package/src/readfos.py
index 80e07f2..ff7885f 100755
--- a/package/src/readfos.py
+++ b/package/src/readfos.py
@@ -271,18 +271,9 @@ class specpol(object):
         """
         data_dump = np.loadtxt(join_path(data_dir, filename), skiprows=1).T
         self.zero(data_dump.shape[1])
-        (
-            self.wav,
-            self.wav_err[:, 0],
-            self.I,
-            self.IQUV_cov[0, 0],
-            self.Q,
-            self.IQUV_cov[1, 1],
-            self.U,
-            self.IQUV_cov[2, 2],
-            self.V,
-            self.IQUV_cov[3, 3],
-        ) = data_dump[:10]
+        (self.wav, self.wav_err[:, 0], self.I, self.IQUV_cov[0, 0], self.Q, self.IQUV_cov[1, 1], self.U, self.IQUV_cov[2, 2], self.V, self.IQUV_cov[3, 3]) = (
+            data_dump[:10]
+        )
         self.wav_err[:, 1] = deepcopy(self.wav_err[:, 0])
         self.bin_edges[:-1], self.bin_edges[-1] = deepcopy(self.wav - self.wav_err[:, 0]), deepcopy(self.wav[-1] + self.wav_err[-1, 1])
         for i in range(4):
@@ -314,11 +305,7 @@ class specpol(object):
                 self.PA_err,
             ]
         ).T
-        np.savetxt(
-            join_path(output_dir, filename + ".txt"),
-            data_dump,
-            header=header,
-        )
+        np.savetxt(join_path(output_dir, filename + ".txt"), data_dump, header=header)
         return join_path(output_dir, filename)
 
     def plot(self, fig=None, ax=None, savename=None, plots_folder=""):