modify shebang

package/FOC_reduction.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python3
 # -*- coding:utf-8 -*-
 """
 Main script where are progressively added the steps for the FOC pipeline reduction.
@@ -30,7 +30,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
         # from lib.deconvolve import from_file_psf
         psf = "gaussian"  # Can be user-defined as well
         # psf = from_file_psf(data_folder+psf_file)
-        psf_FWHM = 3.1
+        psf_FWHM = 1.55
         psf_scale = "px"
         psf_shape = None  # (151, 151)
         iterations = 1
@@ -42,7 +42,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
     # Background estimation
     error_sub_type = "freedman-diaconis"  # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51))
     subtract_error = 1.0
-    display_bkg = True
+    display_bkg = False
 
     # Data binning
     pxsize = 0.05
@@ -51,7 +51,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
 
     # Alignement
     align_center = "center"  # If None will not align the images
-    display_align = True
+    display_align = False
     display_data = False
 
     # Transmittance correction
@@ -66,10 +66,10 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
     rotate_North = True
 
     #  Polarization map output
-    P_cut = 0.99  # if >=1.0 cut on the signal-to-noise else cut on the confidence level in Q, U
+    P_cut = 0.999  # if >=1.0 cut on the signal-to-noise else cut on the confidence level in Q, U
-    SNRi_cut = 1.0  # I measurments with SNR>30, which implies an uncertainty in P of 4.7%.
+    SNRi_cut = 3.0  # I measurments with SNR>30, which implies an uncertainty in P of 4.7%.
     flux_lim = None  # lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None
-    scale_vec = 3
+    scale_vec = 2
     step_vec = 1  # plot all vectors in the array. if step_vec = 2, then every other vector will be plotted if step_vec = 0 then all vectors are displayed at full length
 
     # Pipeline start

package/lib/deconvolve.py

@@ -286,11 +286,13 @@ def richardson_lucy(image, psf, iterations=20, clip=True, filter_epsilon=None):
     image = image.astype(float_type, copy=False)
     psf = psf.astype(float_type, copy=False)
     psf /= psf.sum()
-    im_deconv = image.copy()
+    im_deconv = np.full(image.shape, 0.5, dtype=float_type)
     psf_mirror = np.flip(psf)
 
+    eps = 1e-20
+
     for _ in range(iterations):
-        conv = convolve(im_deconv, psf, mode="same")
+        conv = convolve(im_deconv, psf, mode="same") + eps
         if filter_epsilon:
             relative_blur = np.where(conv < filter_epsilon, 0, image / conv)
         else:

package/lib/plots.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python3
 """
 Library functions for displaying  informations using matplotlib
 
@@ -270,6 +270,7 @@ def polarization_map(
     scale_vec : float, optional
         Pixel length of displayed 100% polarization vector.
         If scale_vec = 2, a vector of 50% polarization will be 1 pixel wide.
+        If scale_vec = 0, all polarization vectors will be displayed at full length.
         Defaults to 2.
     savename : str, optional
         Name of the figure the map should be saved to. If None, the map won't
@@ -412,7 +413,7 @@ def polarization_map(
         ax.set_facecolor("white")
         font_color = "black"
 
-    if display.lower() in ["intensity"]:
+    if display.lower() in ["i", "intensity"]:
         # If no display selected, show intensity map
         display = "i"
         if flux_lim is None:
@@ -427,19 +428,22 @@ def polarization_map(
         levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax
         print("Total flux contour levels : ", levelsI)
         ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5)
-    elif display.lower() in ["pol_flux"]:
+    elif display.lower() in ["pf", "pol_flux"]:
         # Display polarization flux
         display = "pf"
         if flux_lim is None:
             if mask.sum() > 0.0:
                 vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux)
+                pfmax = (stkI[mask] * pol[mask] * convert_flux).max()
             else:
                 vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux)
+                pfmax = (stkI[stkI > 0.0] * pol[stkI > 0.0] * convert_flux).max()
         else:
             vmin, vmax = flux_lim
         im = ax.imshow(stkI * convert_flux * pol, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0)
         fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
-        levelsPf = np.linspace(vmax * 0.01, vmax * 0.99, 10)
+        # levelsPf = np.linspace(0.0175, 0.50, 5) * pfmax
+        levelsPf = np.array([1.73, 13.0, 33.0, 66.0]) / 100.0 * pfmax
         print("Polarized flux contour levels : ", levelsPf)
         ax.contour(stkI * convert_flux * pol, levels=levelsPf, colors="grey", linewidths=0.5)
     elif display.lower() in ["p", "pol", "pol_deg"]:
@@ -548,8 +552,8 @@ def polarization_map(
         arrow_props={"ec": "k", "fc": font_color, "alpha": 1, "lw": 1},
     )
 
-    if display.lower() in ["i", "s_i", "snri", "pf", "p", "pa", "s_p", "snrp", "confp"]:
+    if display.lower() in ["i", "s_i", "snri", "pf", "p", "pa", "s_p", "snrp", "confp"] and step_vec != 0:
-        if step_vec == 0:
+        if scale_vec == -1:
             poldata[np.isfinite(poldata)] = 1.0 / 2.0
             step_vec = 1
             scale_vec = 2.0
@@ -1359,7 +1363,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=None, P_cut=0.99, SNRi_cut=1.0, scale_vec=2.0, step_vec=1, savename=None, **kwargs):
+    def overplot(self, levels=None, 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
@@ -1402,13 +1406,6 @@ class overplot_pol(align_maps):
 
         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} observation from {1:s} overplotted with polarization vectors and Stokes I contours from {2:s}".format(
-                obj, self.other_observer, self.map_observer
-            ),
-            wrap=True,
-        )
-
         # 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
         for key, value in [
@@ -1459,59 +1456,85 @@ class overplot_pol(align_maps):
         )
 
         # Display full size polarization vectors
-        if scale_vec is None:
+        vecstr = ""
-            self.scale_vec = 2.0 * self.px_scale
+        if step_vec != 0:
-            pol[np.isfinite(pol)] = 1.0 / 2.0
+            vecstr = "polarization vectors "
-        else:
+            if scale_vec is None:
-            self.scale_vec = scale_vec * self.px_scale
+                self.scale_vec = 2.0 * self.px_scale
-        self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
+                pol[np.isfinite(pol)] = 1.0 / 2.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)
+            else:
-        self.Q = self.ax_overplot.quiver(
+                self.scale_vec = scale_vec * self.px_scale
-            self.X[::step_vec, ::step_vec],
+            self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
-            self.Y[::step_vec, ::step_vec],
+            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[::step_vec, ::step_vec],
+            self.Q = self.ax_overplot.quiver(
-            self.V[::step_vec, ::step_vec],
+                self.X[::step_vec, ::step_vec],
-            units="xy",
+                self.Y[::step_vec, ::step_vec],
-            angles="uv",
+                self.U[::step_vec, ::step_vec],
-            scale=1.0 / self.scale_vec,
+                self.V[::step_vec, ::step_vec],
-            scale_units="xy",
+                units="xy",
-            pivot="mid",
+                angles="uv",
-            headwidth=0.0,
+                scale=1.0 / self.scale_vec,
-            headlength=0.0,
+                scale_units="xy",
-            headaxislength=0.0,
+                pivot="mid",
-            width=kwargs["width"],
+                headwidth=0.0,
-            linewidth=kwargs["linewidth"],
+                headlength=0.0,
-            color="white",
+                headaxislength=0.0,
-            edgecolor="black",
+                width=kwargs["width"],
-            transform=self.ax_overplot.get_transform(self.wcs_UV),
+                linewidth=kwargs["linewidth"],
-            label="{0:s} polarization map".format(self.map_observer),
+                color="white",
-        )
+                edgecolor="black",
+                transform=self.ax_overplot.get_transform(self.wcs_UV),
+                label="{0:s} polarization map".format(self.map_observer),
+            )
 
-        # Display Stokes I as contours
+        # Display Stokes as contours
-        if levels is None:
+        disptypestr = ""
-            levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(stkI[stkI > 0.0]) * self.map_convert
+        if disptype.lower() == "p":
-        cont_stkI = self.ax_overplot.contour(
+            disptypestr = "polarization degree"
-            stkI * self.map_convert, levels=levels, colors="grey", alpha=0.75, transform=self.ax_overplot.get_transform(self.wcs_UV)
+            if levels is None:
-        )
+                levels = np.array([2.0, 5.0, 10.0, 20.0, 90.0]) / 100.0 * np.max(pol[stkI > 0.0])
-        # self.ax_overplot.clabel(cont_stkI, inline=True, fontsize=5)
+            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)
+            )
+        if disptype.lower() == "pf":
+            disptypestr = "polarized flux"
+            if levels is None:
+                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)
+            )
+        if disptype.lower() == "snri":
+            disptypestr = "Stokes I signal-to-noise"
+            if levels is None:
+                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))
+        else:  # default to intensity contours
+            disptypestr = "Stokes I"
+            if levels is None:
+                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)
+            )
+        # self.ax_overplot.clabel(cont_stk, inline=False, colors="k", fontsize=7)
 
         # Display pixel scale and North direction
-        fontprops = fm.FontProperties(size=16)
+        if step_vec != 0:
-        px_size = self.other_wcs.wcs.get_cdelt()[0] * 3600.0
+            fontprops = fm.FontProperties(size=16)
-        px_sc = AnchoredSizeBar(
+            px_size = self.other_wcs.wcs.get_cdelt()[0] * 3600.0
-            self.ax_overplot.transData,
+            px_sc = AnchoredSizeBar(
-            1.0 / px_size,
+                self.ax_overplot.transData,
-            "1 arcsec",
+                1.0 / px_size,
-            3,
+                "1 arcsec",
-            pad=0.5,
+                3,
-            sep=5,
+                pad=0.5,
-            borderpad=0.5,
+                sep=5,
-            frameon=False,
+                borderpad=0.5,
-            size_vertical=0.005,
+                frameon=False,
-            color=font_color,
+                size_vertical=0.005,
-            fontproperties=fontprops,
+                color=font_color,
-        )
+                fontproperties=fontprops,
-        self.ax_overplot.add_artist(px_sc)
+            )
+            self.ax_overplot.add_artist(px_sc)
         north_dir = AnchoredDirectionArrows(
             self.ax_overplot.transAxes,
             "E",
@@ -1528,20 +1551,21 @@ class overplot_pol(align_maps):
             arrow_props={"ec": "k", "fc": "w", "alpha": 1, "lw": 0.5},
         )
         self.ax_overplot.add_artist(north_dir)
-        pol_sc = AnchoredSizeBar(
+        if step_vec != 0:
-            self.ax_overplot.transData,
+            pol_sc = AnchoredSizeBar(
-            self.scale_vec,
+                self.ax_overplot.transData,
-            r"$P$= 100%",
+                self.scale_vec,
-            4,
+                r"$P$= 100%",
-            pad=0.5,
+                4,
-            sep=5,
+                pad=0.5,
-            borderpad=0.5,
+                sep=5,
-            frameon=False,
+                borderpad=0.5,
-            size_vertical=0.005,
+                frameon=False,
-            color=font_color,
+                size_vertical=0.005,
-            fontproperties=fontprops,
+                color=font_color,
-        )
+                fontproperties=fontprops,
-        self.ax_overplot.add_artist(pol_sc)
+            )
+            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_other,) = self.ax_overplot.plot(*(self.other_wcs.celestial.wcs.crpix - (1.0, 1.0)), "g+")
@@ -1554,15 +1578,21 @@ class overplot_pol(align_maps):
             self.legend_title = r"{0:s} image".format(self.other_observer)
 
         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(
+        if step_vec != 0:
-            (0, 0), (0, 1), arrowstyle="-", fc="w", ec="k", lw=2
+            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
-        labels.append("{0:s} Stokes I contour".format(self.map_observer))
+            )
-        handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=cont_stkI.get_edgecolor()[0]))
+        labels.append("{0:s} {1:s} contour".format(self.map_observer, disptypestr))
+        handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=cont_stk.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
         )
 
+        self.fig_overplot.suptitle(
+            "{0:s} observation from {1:s} overplotted with {2:s} contours from {3:s}".format(obj, self.other_observer, vecstr + disptypestr, self.map_observer),
+            wrap=True,
+        )
+
         if savename is not None:
             if savename[-4:] not in [".png", ".jpg", ".pdf"]:
                 savename += ".pdf"
@@ -1570,10 +1600,10 @@ class overplot_pol(align_maps):
 
         self.fig_overplot.canvas.draw()
 
-    def plot(self, levels=None, P_cut=0.99, SNRi_cut=1.0, scale_vec=2.0, 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, scale_vec=scale_vec, 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):
@@ -3360,8 +3390,17 @@ class pol_map(object):
                 + r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err * 10.0) / 10.0)
                 + str_conf
             )
-            self.str_cut = ""
+            # self.str_cut = ""
-            # self.str_cut = "\n"+r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav, sci_not(I_cut*self.map_convert, I_cut_err*self.map_convert, 2))+"\n"+r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut*100., np.ceil(P_cut_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err*10.)/10.)
+            self.str_cut = (
+                "\n"
+                + r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(
+                    self.pivot_wav, sci_not(I_cut * self.map_convert, I_cut_err * self.map_convert, 2)
+                )
+                + "\n"
+                + r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut * 100.0, np.ceil(P_cut_err * 1000.0) / 10.0)
+                + "\n"
+                + r"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err * 10.0) / 10.0)
+            )
             self.an_int = ax.annotate(
                 self.str_int + self.str_cut,
                 color="white",
@@ -3387,8 +3426,17 @@ class pol_map(object):
                 + r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err * 10.0) / 10.0)
                 + str_conf
             )
-            str_cut = ""
+            # str_cut = ""
-            # str_cut = "\n"+r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav, sci_not(I_cut*self.map_convert, I_cut_err*self.map_convert, 2))+"\n"+r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut*100., np.ceil(P_cut_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err*10.)/10.)
+            str_cut = (
+                "\n"
+                + r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(
+                    self.pivot_wav, sci_not(I_cut * self.map_convert, I_cut_err * self.map_convert, 2)
+                )
+                + "\n"
+                + r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut * 100.0, np.ceil(P_cut_err * 1000.0) / 10.0)
+                + "\n"
+                + r"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err * 10.0) / 10.0)
+            )
             ax.annotate(
                 str_int + str_cut,
                 color="white",

package/lib/query.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python3
+#!/usr/bin/env python3
 # -*- coding:utf-8 -*-
 """
 Library function to query and download datatsets from MAST api.

package/lib/reduction.py

@@ -416,12 +416,12 @@ def deconvolve_array(data_array, headers, psf="gaussian", FWHM=1.0, scale="px",
         FWHM /= pxsize[0].min()
 
     # Define Point-Spread-Function kernel
-    if psf.lower() in ["gauss", "gaussian"]:
+    if isinstance(psf, np.ndarray) and (len(psf.shape) == 2):
+        kernel = psf
+    elif psf.lower() in ["gauss", "gaussian"]:
         if shape is None:
             shape = np.min(data_array[0].shape) - 2, np.min(data_array[0].shape) - 2
         kernel = gaussian_psf(FWHM=FWHM, shape=shape)
-    elif isinstance(psf, np.ndarray) and (len(psf.shape) == 2):
-        kernel = psf
     else:
         raise ValueError("{} is not a valid value for 'psf'".format(psf))
 

package/src/Combine.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python3
 # -*- coding:utf-8 -*-
 from pathlib import Path
 from sys import path as syspath

package/src/emission_center.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python3
 # -*- coding:utf-8 -*-
 from pathlib import Path
 from sys import path as syspath
@@ -41,18 +41,20 @@ def main(infile, target=None, output_dir=None):
         target = Stokes[0].header["TARGNAME"]
 
     fig = figure(figsize=(8, 8.5), layout="constrained")
-    fig, ax = polarization_map(Stokes, P_cut=0.99, step_vec=1, scale_vec=3, display="i", fig=fig, width=0.33, linewidth=0.5)
+    fig, ax = polarization_map(Stokes, P_cut=0.99, step_vec=1, scale_vec=3, display="pf", fig=fig, width=0.33, linewidth=0.5)
 
-    ax.plot(*Stokescenter, marker="+", color="gray", label="Best confidence for center: {0}".format(Stokespos.to_string("hmsdms")))
+    ax.plot(*Stokescenter, marker="+", color="k", lw=3)
+    ax.plot(*Stokescenter, marker="+", color="gray", 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="gray")
-    confcont = ax.contour(Stokesconf, [0.9905], colors="r")
+    # confcont = ax.contour(Stokesconf, [0.9905], colors="r")
     # snr3cont = ax.contour(Stokessnr, [3.0], colors="b", linestyles="dashed")
     # snr4cont = ax.contour(Stokessnr, [4.0], colors="b")
     handles, labels = ax.get_legend_handles_labels()
     labels.append(r"Center $Conf_{99\%}$ contour")
     handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=confcentcont.get_edgecolor()[0]))
-    labels.append(r"Polarization $Conf_{99\%}$ contour")
+    # labels.append(r"Polarization $Conf_{99\%}$ contour")
-    handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=confcont.get_edgecolor()[0]))
+    # handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=confcont.get_edgecolor()[0]))
     # labels.append(r"$SNR_P \geq$ 3  contour")
     # handles.append(Rectangle((0, 0), 1, 1, fill=False, ls="--", ec=snr3cont.get_edgecolor()[0]))
     # labels.append(r"$SNR_P \geq$ 4  contour")

package/src/get_cdelt.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python3
 # -*- coding:utf-8 -*-
 from pathlib import Path
 from sys import path as syspath
@@ -18,7 +18,7 @@ def main(infiles=None):
     from numpy.linalg import eig
 
     if infiles is None:
-        print('Usage: "python get_cdelt.py -f infiles"')
+        print('Usage: "env python3 get_cdelt.py -f infiles"')
         return 1
     prod = [["/".join(filepath.split("/")[:-1]), filepath.split("/")[-1]] for filepath in infiles]
     data_folder = prod[0][0]

package/src/overplot_IC5063.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python3
+#!/usr/bin/env python3
 # -*- coding:utf-8 -*-
 from pathlib import Path
 from sys import path as syspath

package/src/overplot_MRK463E.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python3
+#!/usr/bin/env python3
 # -*- coding:utf-8 -*-
 from pathlib import Path
 from sys import path as syspath
@@ -26,18 +26,21 @@ Radio = fits.open("./data/MRK463E/EMERLIN/Voorwerpjes_1356+1822_1356+1822_unifor
 # B.plot(levels=levels, P_cut=0.99, SNRi_cut=1.0, scale_vec=5, norm=LogNorm(8.5e-18, 2.5e-15), savename="./plots/MRK463E/IR_overplot.pdf")
 # B.write_to(path1="./data/MRK463E/FOC_data_WFPC.fits", path2="./data/MRK463E/WFPC_data.fits", suffix="aligned")
 
-levels = np.array([0.8, 2, 5, 10, 20, 50]) / 100.0 * Stokes_UV[0].header["photflam"]
+# levels = np.array([0.8, 2, 5, 10, 20, 50]) / 100.0 * Stokes_UV[0].header["photflam"]
+levels = np.array([5, 10, 20, 50])
 C = overplot_pol(Stokes_UV, Radio, norm=LogNorm())
 C.other_im.set(norm=LogNorm(1e-4, 2e-2))
 C.plot(
     levels=levels,
-    P_cut=0.999,
+    P_cut=0.99,
-    SNRi_cut=5.0,
+    SNRi_cut=1.0,
+    step_vec=0,
     scale_vec=3,
     norm=LogNorm(1e-4, 2e-2),
     cmap="inferno_r",
     width=0.5,
     linewidth=0.5,
-    savename="./plots/MRK463E/EMERLIN_overplot.pdf",
+    disptype="snri",
+    savename="./plots/MRK463E/EMERLIN_snri_overplot.pdf",
 )
 C.write_to(path1="./data/MRK463E/FOC_data_EMERLIN.fits", path2="./data/MRK463E/EMERLIN_data.fits", suffix="aligned")

package/src/overplot_NGC1068.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python3
+#!/usr/bin/env python3
 # -*- coding:utf-8 -*-
 from pathlib import Path
 from sys import path as syspath