fix pol_map.set_data_mask and some keyword saving

finish rebase display on main

README.md

@@ -1,2 +1,6 @@
 # FOC_Reduction
 FOC reduction pipeline
+
+TODO:
+    - Add all polarimetry capables instruments from HST (starting with FOS and ACS)
+    - Build science case for future UV polarimeters (AGN outflow geometry, dust scattering, torus outline ?)

license.md

@@ -0,0 +1,10 @@
+MIT License
+
+Copyright (c) 2024 Thibault Barnouin
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

package/FOC_reduction.py

@@ -1,9 +1,14 @@
-#!/usr/bin/python
+#!/usr/bin/env python3
 # -*- coding:utf-8 -*-
 """
 Main script where are progressively added the steps for the FOC pipeline reduction.
 """
 
+from pathlib import Path
+from sys import path as syspath
+
+syspath.append(str(Path(__file__).parent.parent))
+
 # Project libraries
 from copy import deepcopy
 from os import system
@@ -15,6 +20,7 @@ import lib.reduction as proj_red  # Functions used in reduction pipeline
 import numpy as np
 from lib.utils import princ_angle, sci_not
 from matplotlib.colors import LogNorm
+from astropy.wcs import WCS
 
 
 def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=False, interactive=False):
@@ -25,7 +31,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
@@ -35,9 +41,9 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
     display_crop = False
 
     # Background estimation
-    error_sub_type = "freedman-diaconis"  # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51))
+    error_sub_type = "scott"  # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51))
-    subtract_error = 1.0
+    subtract_error = 1.50
-    display_bkg = False
+    display_bkg = True
 
     # Data binning
     pxsize = 0.05
@@ -54,29 +60,20 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
 
     # Smoothing
     smoothing_function = "combine"  # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine
-    smoothing_FWHM = 0.10  # If None, no smoothing is done
+    smoothing_FWHM = 0.075  # If None, no smoothing is done
     smoothing_scale = "arcsec"  # pixel or arcsec
 
     # Rotation
+    rotate_data = False
     rotate_North = True
 
     #  Polarization map output
-    SNRp_cut = 3.0  # P measurments with SNR>3
+    P_cut = 3  # 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%.
     flux_lim = None  # lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None
     scale_vec = 5
     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
 
-    # Adaptive binning
-    # in order to perfrom optimal binning, there are several steps to follow:
-    # 1. Load the data again and preserve the full images
-    # 2. Skip the cropping step but use the same error and background estimation
-    # 3. Use the same alignment as the routine
-    # 4. Skip the rebinning step
-    # 5. Calulate the Stokes parameters without smoothing
-    optimal_binning = False
-    optimize = False
-
     # Pipeline start
 
     # Step 1:
@@ -121,278 +118,6 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
     if align_center is None:
         figtype = "_".join([figtype, "not_aligned"] if figtype != "" else ["not_aligned"])
 
-    if optimal_binning:
-        from lib.background import subtract_bkg
-
-        options = {"optimize": optimize, "optimal_binning": True}
-
-        # Step 1: Load the data again and preserve the full images
-        _data_array, _headers = deepcopy(data_array), deepcopy(headers)  # Preserve full images
-        _data_mask = np.ones(_data_array[0].shape, dtype=bool)
-
-        # Step 2: Skip the cropping step but use the same error and background estimation (I don't understand why this is wrong)
-        data_array, error_array, headers = proj_red.crop_array(
-            data_array, headers, step=5, null_val=0.0, inside=True, display=display_crop, savename=figname, plots_folder=plots_folder
-        )
-        data_mask = np.ones(data_array[0].shape, dtype=bool)
-
-        background = None
-        _, _, _, background, error_bkg = proj_red.get_error(
-            data_array,
-            headers,
-            error_array,
-            data_mask=data_mask,
-            sub_type=error_sub_type,
-            subtract_error=subtract_error,
-            display=display_bkg,
-            savename="_".join([figname, "errors"]),
-            plots_folder=plots_folder,
-            return_background=True,
-        )
-
-        # _background is the same as background, but for the optimal binning
-        _background = None
-        _data_array, _error_array, _ = proj_red.get_error(
-            _data_array,
-            _headers,
-            error_array=None,
-            data_mask=_data_mask,
-            sub_type=error_sub_type,
-            subtract_error=False,
-            display=display_bkg,
-            savename="_".join([figname, "errors"]),
-            plots_folder=plots_folder,
-            return_background=False,
-        )
-        _error_bkg = np.ones_like(_data_array) * error_bkg[:, 0, 0, np.newaxis, np.newaxis]
-        _data_array, _error_array, _background, _ = subtract_bkg(_data_array, _error_array, _data_mask, background, _error_bkg)
-
-        # Step 3: Align and rescale images with oversampling. (has to disable croping in align_data function)
-        _data_array, _error_array, _headers, _, shifts, error_shifts = proj_red.align_data(
-            _data_array,
-            _headers,
-            error_array=_error_array,
-            background=_background,
-            upsample_factor=10,
-            ref_center=align_center,
-            return_shifts=True,
-            optimal_binning=True,
-        )
-        print("Image shifts: {} \nShifts uncertainty: {}".format(shifts, error_shifts))
-        _data_mask = np.ones(_data_array[0].shape, dtype=bool)
-
-        # Step 4: Compute Stokes I, Q, U
-        _background = np.array([np.array(bkg).reshape(1, 1) for bkg in _background])
-        _background_error = np.array(
-            [
-                np.array(
-                    np.sqrt(
-                        (bkg - _background[np.array([h["filtnam1"] == head["filtnam1"] for h in _headers], dtype=bool)].mean()) ** 2
-                        / np.sum([h["filtnam1"] == head["filtnam1"] for h in _headers])
-                    )
-                ).reshape(1, 1)
-                for bkg, head in zip(_background, _headers)
-            ]
-        )
-
-        _I_stokes, _Q_stokes, _U_stokes, _Stokes_cov, _header_stokes = proj_red.compute_Stokes(
-            _data_array, _error_array, _data_mask, _headers, FWHM=None, scale=smoothing_scale, smoothing=smoothing_function, transmitcorr=transmitcorr
-        )
-        _I_bkg, _Q_bkg, _U_bkg, _S_cov_bkg, _header_bkg = proj_red.compute_Stokes(
-            _background,
-            _background_error,
-            np.array(True).reshape(1, 1),
-            _headers,
-            FWHM=None,
-            scale=smoothing_scale,
-            smoothing=smoothing_function,
-            transmitcorr=False,
-        )
-
-        # Step 5: Compute polarimetric parameters (polarization degree and angle).
-        _P, _debiased_P, _s_P, _s_P_P, _PA, _s_PA, _s_PA_P = proj_red.compute_pol(_I_stokes, _Q_stokes, _U_stokes, _Stokes_cov, _header_stokes)
-        _P_bkg, _debiased_P_bkg, _s_P_bkg, _s_P_P_bkg, _PA_bkg, _s_PA_bkg, _s_PA_P_bkg = proj_red.compute_pol(_I_bkg, _Q_bkg, _U_bkg, _S_cov_bkg, _header_bkg)
-
-        # Step 6: Save image to FITS.
-        figname = "_".join([figname, figtype]) if figtype != "" else figname
-        _Stokes_hdul = proj_fits.save_Stokes(
-            _I_stokes,
-            _Q_stokes,
-            _U_stokes,
-            _Stokes_cov,
-            _P,
-            _debiased_P,
-            _s_P,
-            _s_P_P,
-            _PA,
-            _s_PA,
-            _s_PA_P,
-            _header_stokes,
-            _data_mask,
-            figname,
-            data_folder=data_folder,
-            return_hdul=True,
-        )
-
-        # Step 6:
-        _data_mask = _Stokes_hdul["data_mask"].data.astype(bool)
-        print(
-            "F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(
-                _header_stokes["PHOTPLAM"],
-                *sci_not(
-                    _Stokes_hdul[0].data[_data_mask].sum() * _header_stokes["PHOTFLAM"],
-                    np.sqrt(_Stokes_hdul[3].data[0, 0][_data_mask].sum()) * _header_stokes["PHOTFLAM"],
-                    2,
-                    out=int,
-                ),
-            )
-        )
-        print("P_int = {0:.1f} ± {1:.1f} %".format(_header_stokes["p_int"] * 100.0, np.ceil(_header_stokes["sP_int"] * 1000.0) / 10.0))
-        print("PA_int = {0:.1f} ± {1:.1f} °".format(princ_angle(_header_stokes["pa_int"]), princ_angle(np.ceil(_header_stokes["sPA_int"] * 10.0) / 10.0)))
-        #  Background values
-        print(
-            "F_bkg({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(
-                _header_stokes["PHOTFLAM"],
-                *sci_not(_I_bkg[0, 0] * _header_stokes["PHOTFLAM"], np.sqrt(_S_cov_bkg[0, 0][0, 0]) * _header_stokes["PHOTFLAM"], 2, out=int),
-            )
-        )
-        print("P_bkg = {0:.1f} ± {1:.1f} %".format(_debiased_P_bkg[0, 0] * 100.0, np.ceil(_s_P_bkg[0, 0] * 1000.0) / 10.0))
-        print("PA_bkg = {0:.1f} ± {1:.1f} °".format(princ_angle(_PA_bkg[0, 0]), princ_angle(np.ceil(_s_PA_bkg[0, 0] * 10.0) / 10.0)))
-
-        #  Plot polarization map (Background is either total Flux, Polarization degree or Polarization degree error).
-        if pxscale.lower() not in ["full", "integrate"] and not interactive:
-            proj_plots.polarization_map(
-                deepcopy(_Stokes_hdul),
-                _data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                vec_scale=scale_vec,
-                savename="_".join([figname]),
-                plots_folder=plots_folder,
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(_Stokes_hdul),
-                _data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                vec_scale=scale_vec,
-                savename="_".join([figname, "I"]),
-                plots_folder=plots_folder,
-                display="Intensity",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(_Stokes_hdul),
-                _data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                vec_scale=scale_vec,
-                savename="_".join([figname, "P_flux"]),
-                plots_folder=plots_folder,
-                display="Pol_Flux",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(_Stokes_hdul),
-                _data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                vec_scale=scale_vec,
-                savename="_".join([figname, "P"]),
-                plots_folder=plots_folder,
-                display="Pol_deg",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(_Stokes_hdul),
-                _data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                vec_scale=scale_vec,
-                savename="_".join([figname, "PA"]),
-                plots_folder=plots_folder,
-                display="Pol_ang",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(_Stokes_hdul),
-                _data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                vec_scale=scale_vec,
-                savename="_".join([figname, "I_err"]),
-                plots_folder=plots_folder,
-                display="I_err",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(_Stokes_hdul),
-                _data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                vec_scale=scale_vec,
-                savename="_".join([figname, "P_err"]),
-                plots_folder=plots_folder,
-                display="Pol_deg_err",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(_Stokes_hdul),
-                _data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                vec_scale=scale_vec,
-                savename="_".join([figname, "SNRi"]),
-                plots_folder=plots_folder,
-                display="SNRi",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(_Stokes_hdul),
-                _data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                vec_scale=scale_vec,
-                savename="_".join([figname, "SNRp"]),
-                plots_folder=plots_folder,
-                display="SNRp",
-                **options,
-            )
-        elif not interactive:
-            proj_plots.polarization_map(
-                deepcopy(_Stokes_hdul),
-                _data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                savename=figname,
-                plots_folder=plots_folder,
-                display="integrate",
-                **options,
-            )
-        elif pxscale.lower() not in ["full", "integrate"]:
-            proj_plots.pol_map(_Stokes_hdul, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim)
-
-    else:
-        options = {"optimize": optimize, "optimal_binning": False}
     #  Crop data to remove outside blank margins.
     data_array, error_array, headers = proj_red.crop_array(
         data_array, headers, step=5, null_val=0.0, inside=True, display=display_crop, savename=figname, plots_folder=plots_folder
@@ -401,13 +126,11 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
 
     #  Deconvolve data using Richardson-Lucy iterative algorithm with a gaussian PSF of given FWHM.
     if deconvolve:
-            data_array = proj_red.deconvolve_array(
+        data_array = proj_red.deconvolve_array(data_array, headers, psf=psf, FWHM=psf_FWHM, scale=psf_scale, shape=psf_shape, iterations=iterations, algo=algo)
-                data_array, headers, psf=psf, FWHM=psf_FWHM, scale=psf_scale, shape=psf_shape, iterations=iterations, algo=algo
-            )
 
     #  Estimate error from data background, estimated from sub-image of desired sub_shape.
     background = None
-        data_array, error_array, headers, background, error_bkg = proj_red.get_error(
+    data_array, error_array, headers, background = proj_red.get_error(
         data_array,
         headers,
         error_array,
@@ -420,31 +143,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
         return_background=True,
     )
 
-        #  Align and rescale images with oversampling.
-        data_array, error_array, headers, data_mask, shifts, error_shifts = proj_red.align_data(
-            data_array, headers, error_array=error_array, background=background, upsample_factor=10, ref_center=align_center, return_shifts=True
-        )
-
-        if display_align:
-            print("Image shifts: {} \nShifts uncertainty: {}".format(shifts, error_shifts))
-            proj_plots.plot_obs(
-                data_array,
-                headers,
-                savename="_".join([figname, str(align_center)]),
-                plots_folder=plots_folder,
-                norm=LogNorm(
-                    vmin=data_array[data_array > 0.0].min() * headers[0]["photflam"], vmax=data_array[data_array > 0.0].max() * headers[0]["photflam"]
-                ),
-            )
-
-        #  Rebin data to desired pixel size.
-        if (pxsize is not None) and not (pxsize == 1 and pxscale.lower() in ["px", "pixel", "pixels"]):
-            data_array, error_array, headers, Dxy, data_mask = proj_red.rebin_array(
-                data_array, error_array, headers, pxsize=pxsize, scale=pxscale, operation=rebin_operation, data_mask=data_mask
-            )
-
     # Rotate data to have same orientation
-        rotate_data = np.unique([np.round(float(head["ORIENTAT"]), 3) for head in headers]).size != 1
+    rotate_data = rotate_data or np.unique([np.round(float(head["ORIENTAT"]), 3) for head in headers]).size != 1
     if rotate_data:
         ang = np.mean([head["ORIENTAT"] for head in headers])
         for head in headers:
@@ -460,6 +160,42 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
                     vmin=data_array[data_array > 0.0].min() * headers[0]["photflam"], vmax=data_array[data_array > 0.0].max() * headers[0]["photflam"]
                 ),
             )
+    #  Align and rescale images with oversampling.
+    data_array, error_array, headers, data_mask, shifts, error_shifts = proj_red.align_data(
+        data_array,
+        headers,
+        error_array=error_array,
+        data_mask=data_mask,
+        background=background,
+        upsample_factor=10,
+        ref_center=align_center,
+        return_shifts=True,
+    )
+
+    if display_align:
+        print("Image shifts: {} \nShifts uncertainty: {}".format(shifts, error_shifts))
+        proj_plots.plot_obs(
+            data_array,
+            headers,
+            shifts=shifts,
+            savename="_".join([figname, str(align_center)]),
+            plots_folder=plots_folder,
+            norm=LogNorm(vmin=data_array[data_array > 0.0].min() * headers[0]["photflam"], vmax=data_array[data_array > 0.0].max() * headers[0]["photflam"]),
+        )
+
+    flux_data, flux_error, flux_mask, flux_head = (
+        deepcopy(data_array.sum(axis=0)),
+        deepcopy(np.sqrt(np.sum(error_array**2, axis=0))),
+        deepcopy(data_mask),
+        deepcopy(headers[0]),
+    )
+    flux_head["EXPTIME"] = np.sum([head["EXPTIME"] for head in headers])
+
+    #  Rebin data to desired pixel size.
+    if (pxsize is not None) and not (pxsize == 1 and pxscale.lower() in ["px", "pixel", "pixels"]):
+        data_array, error_array, headers, Dxy, data_mask = proj_red.rebin_array(
+            data_array, error_array, headers, pxsize=pxsize, scale=pxscale, operation=rebin_operation, data_mask=data_mask
+        )
 
     # Plot array for checking output
     if display_data and pxscale.lower() not in ["full", "integrate"]:
@@ -468,9 +204,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
             headers,
             savename="_".join([figname, "rebin"]),
             plots_folder=plots_folder,
-                norm=LogNorm(
+            norm=LogNorm(vmin=data_array[data_array > 0.0].min() * headers[0]["photflam"], vmax=data_array[data_array > 0.0].max() * headers[0]["photflam"]),
-                    vmin=data_array[data_array > 0.0].min() * headers[0]["photflam"], vmax=data_array[data_array > 0.0].max() * headers[0]["photflam"]
-                ),
         )
 
     background = np.array([np.array(bkg).reshape(1, 1) for bkg in background])
@@ -496,14 +230,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
         data_array, error_array, data_mask, headers, FWHM=smoothing_FWHM, scale=smoothing_scale, smoothing=smoothing_function, transmitcorr=transmitcorr
     )
     I_bkg, Q_bkg, U_bkg, S_cov_bkg, header_bkg = proj_red.compute_Stokes(
-            background,
+        background, background_error, np.array(True).reshape(1, 1), headers, FWHM=None, scale=smoothing_scale, smoothing=smoothing_function, transmitcorr=False
-            background_error,
-            np.array(True).reshape(1, 1),
-            headers,
-            FWHM=None,
-            scale=smoothing_scale,
-            smoothing=smoothing_function,
-            transmitcorr=False,
     )
 
     # Step 3:
@@ -515,6 +242,10 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
         I_bkg, Q_bkg, U_bkg, S_cov_bkg, data_mask_bkg, header_bkg = proj_red.rotate_Stokes(
             I_bkg, Q_bkg, U_bkg, S_cov_bkg, np.array(True).reshape(1, 1), header_bkg, SNRi_cut=None
         )
+        flux_data, flux_error, flux_mask, flux_head = proj_red.rotate_data(np.array([flux_data]), np.array([flux_error]), flux_mask, [flux_head])
+        flux_data, flux_error, flux_head = flux_data[0], flux_error[0], flux_head[0]
+    elif not rotate_data:
+        figname += "_noROT"
 
     # Compute polarimetric parameters (polarization degree and angle).
     P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, header_stokes)
@@ -522,7 +253,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
 
     # Step 4:
     # Save image to FITS.
-        figname = "_".join([figname, figtype]) if figtype != "" else figname
+    savename = "_".join([figname, figtype]) if figtype != "" else figname
     Stokes_hdul = proj_fits.save_Stokes(
         I_stokes,
         Q_stokes,
@@ -537,32 +268,29 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
         s_PA_P,
         header_stokes,
         data_mask,
-            figname,
+        savename,
         data_folder=data_folder,
         return_hdul=True,
+        flux_data=np.array([flux_data, flux_error, flux_mask]),
+        flux_head=flux_head,
     )
     outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"]))
 
     # Step 5:
     # crop to desired region of interest (roi)
     if crop:
-            figname += "_crop"
+        savename += "_crop"
         stokescrop = proj_plots.crop_Stokes(deepcopy(Stokes_hdul), norm=LogNorm())
         stokescrop.crop()
         stokescrop.write_to("/".join([data_folder, figname + ".fits"]))
-            Stokes_hdul, header_stokes = stokescrop.hdul_crop, [dataset.header for dataset in stokescrop.hdul_crop]
+        Stokes_hdul, header_stokes = stokescrop.hdul_crop, stokescrop.hdul_crop[0].header
-            outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"]))
+        outfiles.append("/".join([data_folder, Stokes_hdul["I_STOKES"].header["FILENAME"] + ".fits"]))
 
     data_mask = Stokes_hdul["data_mask"].data.astype(bool)
     print(
         "F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(
-                header_stokes["PHOTPLAM"],
+            flux_head["PHOTPLAM"],
-                *sci_not(
+            *sci_not(flux_data[flux_mask].sum() * flux_head["PHOTFLAM"], np.sqrt(np.sum(flux_error[flux_mask] ** 2)) * flux_head["PHOTFLAM"], 2, out=int),
-                    Stokes_hdul[0].data[data_mask].sum() * header_stokes["PHOTFLAM"],
-                    np.sqrt(Stokes_hdul[3].data[0, 0][data_mask].sum()) * header_stokes["PHOTFLAM"],
-                    2,
-                    out=int,
-                ),
         )
     )
     print("P_int = {0:.1f} ± {1:.1f} %".format(header_stokes["p_int"] * 100.0, np.ceil(header_stokes["sP_int"] * 1000.0) / 10.0))
@@ -570,8 +298,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
     #  Background values
     print(
         "F_bkg({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(
-                header_stokes["PHOTPLAM"],
+            header_stokes["photplam"], *sci_not(I_bkg[0, 0] * header_stokes["photflam"], np.sqrt(S_cov_bkg[0, 0][0, 0]) * header_stokes["photflam"], 2, out=int)
-                *sci_not(I_bkg[0, 0] * header_stokes["PHOTPLAM"], np.sqrt(S_cov_bkg[0, 0][0, 0]) * header_stokes["PHOTPLAM"], 2, out=int),
         )
     )
     print("P_bkg = {0:.1f} ± {1:.1f} %".format(debiased_P_bkg[0, 0] * 100.0, np.ceil(s_P_bkg[0, 0] * 1000.0) / 10.0))
@@ -581,132 +308,39 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
         proj_plots.polarization_map(
             deepcopy(Stokes_hdul),
             data_mask,
-                SNRp_cut=SNRp_cut,
+            P_cut=P_cut,
             SNRi_cut=SNRi_cut,
             flux_lim=flux_lim,
             step_vec=step_vec,
             scale_vec=scale_vec,
-                savename="_".join([figname]),
-                plots_folder=plots_folder,
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(Stokes_hdul),
-                data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                scale_vec=scale_vec,
-                savename="_".join([figname, "I"]),
-                plots_folder=plots_folder,
-                display="Intensity",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(Stokes_hdul),
-                data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                scale_vece=scale_vec,
-                savename="_".join([figname, "P_flux"]),
-                plots_folder=plots_folder,
-                display="Pol_Flux",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(Stokes_hdul),
-                data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                scale_vec=scale_vec,
-                savename="_".join([figname, "P"]),
-                plots_folder=plots_folder,
-                display="Pol_deg",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(Stokes_hdul),
-                data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                scale_vec=scale_vec,
-                savename="_".join([figname, "PA"]),
-                plots_folder=plots_folder,
-                display="Pol_ang",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(Stokes_hdul),
-                data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                scale_vec=scale_vec,
-                savename="_".join([figname, "I_err"]),
-                plots_folder=plots_folder,
-                display="I_err",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(Stokes_hdul),
-                data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                scale_vec=scale_vec,
-                savename="_".join([figname, "P_err"]),
-                plots_folder=plots_folder,
-                display="Pol_deg_err",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(Stokes_hdul),
-                data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                scale_vec=scale_vec,
-                savename="_".join([figname, "SNRi"]),
-                plots_folder=plots_folder,
-                display="SNRi",
-                **options,
-            )
-            proj_plots.polarization_map(
-                deepcopy(Stokes_hdul),
-                data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
-                flux_lim=flux_lim,
-                step_vec=step_vec,
-                scale_vec=scale_vec,
-                savename="_".join([figname, "SNRp"]),
-                plots_folder=plots_folder,
-                display="SNRp",
-                **options,
-            )
-        elif not interactive:
-            proj_plots.polarization_map(
-                deepcopy(Stokes_hdul),
-                data_mask,
-                SNRp_cut=SNRp_cut,
-                SNRi_cut=SNRi_cut,
             savename=figname,
             plots_folder=plots_folder,
-                display="integrate",
+        )
-                **options,
+        for figtype, figsuffix in zip(
+            ["FluxDensity", "Intensity", "Pol_flux", "Pol_deg", "Pol_ang", "I_err", "P_err", "SNRi", "SNRp", "confp"],
+            ["F", "I", "P_flux", "P", "PA", "I_err", "P_err", "SNRi", "SNRp", "confP"],
+        ):
+            try:
+                proj_plots.polarization_map(
+                    deepcopy(Stokes_hdul),
+                    data_mask,
+                    P_cut=P_cut,
+                    SNRi_cut=SNRi_cut,
+                    flux_lim=flux_lim,
+                    step_vec=step_vec,
+                    scale_vec=scale_vec,
+                    savename="_".join([savename, figsuffix]),
+                    plots_folder=plots_folder,
+                    display=figtype,
+                )
+            except ValueError:
+                pass
+    elif not interactive:
+        proj_plots.polarization_map(
+            deepcopy(Stokes_hdul), data_mask, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=savename, plots_folder=plots_folder, display="integrate"
         )
     elif pxscale.lower() not in ["full", "integrate"]:
-            proj_plots.pol_map(Stokes_hdul, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim)
+        proj_plots.pol_map(Stokes_hdul, P_cut=P_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, scale_vec=scale_vec, flux_lim=flux_lim)
 
     return outfiles
 

package/__init__.py

@@ -1,3 +1,2 @@
-from . import lib
+from .lib import *
-from . import src
+from .FOC_reduction import main
-from . import FOC_reduction

package/lib/background.py

@@ -18,7 +18,6 @@ import matplotlib.dates as mdates
 import matplotlib.pyplot as plt
 import numpy as np
 from astropy.time import Time
-from lib.plots import plot_obs
 from matplotlib.colors import LogNorm
 from matplotlib.patches import Rectangle
 from scipy.optimize import curve_fit
@@ -85,12 +84,20 @@ def display_bkg(data, background, std_bkg, headers, histograms=None, binning=Non
                 label=headers[i]["filtnam1"] + " (Obs " + str(filt_obs[headers[i]["filtnam1"]]) + ")",
             )
             ax_h.plot([background[i] * convert_flux[i], background[i] * convert_flux[i]], [hist.min(), hist.max()], "x--", color="C{0:d}".format(i), alpha=0.8)
+            if i == 0:
+                xmin, xmax = np.min(np.array(bins)[np.array(hist) > 5e1]) * convert_flux[0], np.max(np.array(bins)[np.array(hist) > 5e1]) * convert_flux[0]
+            else:
+                xmin, xmax = (
+                    min(xmin, np.min(np.array(bins)[np.array(hist) > 5e1]) * convert_flux[0]),
+                    max(xmax, np.max(np.array(bins)[np.array(hist) > 5e1]) * convert_flux[0]),
+                )
             if coeff is not None:
                 # ax_h.plot(bins*convert_flux[i], gausspol(bins, *coeff[i]), '--', color="C{0:d}".format(i), alpha=0.8)
                 ax_h.plot(bins * convert_flux[i], gauss(bins, *coeff[i]), "--", color="C{0:d}".format(i), alpha=0.8)
         ax_h.set_xscale("log")
-        ax_h.set_ylim([0.0, np.max([hist.max() for hist in histograms])])
+        ax_h.set_yscale("log")
-        ax_h.set_xlim([np.min(background * convert_flux) * 1e-2, np.max(background * convert_flux) * 1e2])
+        ax_h.set_ylim([5e1, np.max([hist.max() for hist in histograms])])
+        ax_h.set_xlim([max(xmin, np.min(background * convert_flux) * 1e-2), min(xmax, np.max(background * convert_flux) * 1e2)])
         ax_h.set_xlabel(r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
         ax_h.set_ylabel(r"Number of pixels in bin")
         ax_h.set_title("Histogram for each observation")
@@ -126,7 +133,9 @@ def display_bkg(data, background, std_bkg, headers, histograms=None, binning=Non
     ax2.set(xlabel="pixel offset", ylabel="pixel offset", aspect="equal")
 
     fig2.subplots_adjust(hspace=0, wspace=0, right=1.0)
-    fig2.colorbar(im2, ax=ax2, location="right", aspect=50, pad=0.025, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
+    fig2.colorbar(im2, ax=ax2, location="right", shrink=0.60, aspect=50, pad=0.025, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
+
+    from .plots import plot_obs
 
     if savename is not None:
         this_savename = deepcopy(savename)
@@ -252,20 +261,25 @@ def bkg_fit(data, error, mask, headers, subtract_error=True, display=False, save
         weights /= weights.sum()
 
         bkg = np.sum(weights * (coeff[:, 1] + np.abs(coeff[:, 2]) * subtract_error))
+
         error_bkg[i] *= bkg
 
+        n_error_array[i] = np.sqrt(n_error_array[i] ** 2 + error_bkg[i] ** 2)
+
+        # Substract background
+        if np.abs(subtract_error) > 0:
+            n_data_array[i][mask] = n_data_array[i][mask] - bkg
+            n_data_array[i][np.logical_and(mask, n_data_array[i] <= 1e-3 * bkg)] = 1e-3 * bkg
+
         std_bkg[i] = image[np.abs(image - bkg) / bkg < 1.0].std()
         background[i] = bkg
 
-    if np.abs(subtract_error) > 0:
-        n_data_array, n_error_array, background, error_bkg = subtract_bkg(n_data_array, n_error_array, mask, background, error_bkg)
-
     if display:
         display_bkg(data, background, std_bkg, headers, histograms=histograms, binning=binning, coeff=coeff, savename=savename, plots_folder=plots_folder)
-    return n_data_array, n_error_array, headers, background, error_bkg
+    return n_data_array, n_error_array, headers, background
 
 
-def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, display=False, savename=None, plots_folder=""):
+def bkg_hist(data, error, mask, headers, n_bins=None, subtract_error=True, display=False, savename=None, plots_folder=""):
     """
     ----------
     Inputs:
@@ -281,7 +295,7 @@ def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, dis
     sub_type : str or int, optional
         If str, statistic rule to be used for the number of bins in counts/s.
         If int, number of bins for the counts/s histogram.
-        Defaults to "Freedman-Diaconis".
+        Defaults to "Scott".
     subtract_error : float or bool, optional
         If float, factor to which the estimated background should be multiplied
         If False the background is not subtracted.
@@ -320,29 +334,15 @@ def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, dis
     for i, image in enumerate(data):
         # Compute the Count-rate histogram for the image
         n_mask = np.logical_and(mask, image > 0.0)
-        if sub_type is not None:
-            if isinstance(sub_type, int):
-                n_bins = sub_type
-            elif sub_type.lower() in ["sqrt"]:
-                n_bins = np.fix(np.sqrt(image[n_mask].size)).astype(int)  # Square-root
-            elif sub_type.lower() in ["sturges"]:
-                n_bins = np.ceil(np.log2(image[n_mask].size)).astype(int) + 1  # Sturges
-            elif sub_type.lower() in ["rice"]:
-                n_bins = 2 * np.fix(np.power(image[n_mask].size, 1 / 3)).astype(int)  # Rice
-            elif sub_type.lower() in ["scott"]:
-                n_bins = np.fix((image[n_mask].max() - image[n_mask].min()) / (3.5 * image[n_mask].std() / np.power(image[n_mask].size, 1 / 3))).astype(
-                    int
-                )  # Scott
-            else:
-                n_bins = np.fix(
-                    (image[n_mask].max() - image[n_mask].min())
-                    / (2 * np.subtract(*np.percentile(image[n_mask], [75, 25])) / np.power(image[n_mask].size, 1 / 3))
-                ).astype(int)  # Freedman-Diaconis
-        else:
-            n_bins = np.fix(
-                (image[n_mask].max() - image[n_mask].min()) / (2 * np.subtract(*np.percentile(image[n_mask], [75, 25])) / np.power(image[n_mask].size, 1 / 3))
-            ).astype(int)  # Freedman-Diaconis
 
+        if not isinstance(n_bins, int) and n_bins not in ["auto", "fd", "doane", "scott", "stone", "rice", "sturges", "sqrt"]:
+            match n_bins.lower():
+                case "square-root" | "squareroot":
+                    n_bins = "sqrt"
+                case "freedman-diaconis" | "freedmandiaconis":
+                    n_bins = "fd"
+                case _:
+                    n_bins = "scott"
         hist, bin_edges = np.histogram(np.log(image[n_mask]), bins=n_bins)
         histograms.append(hist)
         binning.append(np.exp(bin_centers(bin_edges)))
@@ -355,19 +355,23 @@ def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, dis
         # popt, pcov = curve_fit(gausspol, binning[-1], hist, p0=p0)
         popt, pcov = curve_fit(gauss, binning[-1], hist, p0=p0)
         coeff.append(popt)
-
         bkg = popt[1] + np.abs(popt[2]) * subtract_error
+
         error_bkg[i] *= bkg
 
+        n_error_array[i] = np.sqrt(n_error_array[i] ** 2 + error_bkg[i] ** 2)
+
+        # Substract background
+        if np.abs(subtract_error) > 0:
+            n_data_array[i][mask] = n_data_array[i][mask] - bkg
+            n_data_array[i][np.logical_and(mask, n_data_array[i] <= 1e-3 * bkg)] = 1e-3 * bkg
+
         std_bkg[i] = image[np.abs(image - bkg) / bkg < 1.0].std()
         background[i] = bkg
 
-    if np.abs(subtract_error) > 0:
-        n_data_array, n_error_array, background, error_bkg = subtract_bkg(n_data_array, n_error_array, mask, background, error_bkg)
-
     if display:
         display_bkg(data, background, std_bkg, headers, histograms=histograms, binning=binning, coeff=coeff, savename=savename, plots_folder=plots_folder)
-    return n_data_array, n_error_array, headers, background, error_bkg
+    return n_data_array, n_error_array, headers, background
 
 
 def bkg_mini(data, error, mask, headers, sub_shape=(15, 15), subtract_error=True, display=False, savename=None, plots_folder=""):
@@ -449,28 +453,19 @@ def bkg_mini(data, error, mask, headers, sub_shape=(15, 15), subtract_error=True
         # Compute error : root mean square of the background
         sub_image = image[minima[0] : minima[0] + sub_shape[0], minima[1] : minima[1] + sub_shape[1]]
         # bkg =  np.std(sub_image)    # Previously computed using standard deviation over the background
-
         bkg = np.sqrt(np.sum(sub_image**2) / sub_image.size) * subtract_error if subtract_error > 0 else np.sqrt(np.sum(sub_image**2) / sub_image.size)
         error_bkg[i] *= bkg
 
+        n_error_array[i] = np.sqrt(n_error_array[i] ** 2 + error_bkg[i] ** 2)
+
+        # Substract background
+        if np.abs(subtract_error) > 0:
+            n_data_array[i][mask] = n_data_array[i][mask] - bkg
+            n_data_array[i][np.logical_and(mask, n_data_array[i] <= 1e-3 * bkg)] = 1e-3 * bkg
+
         std_bkg[i] = image[np.abs(image - bkg) / bkg < 1.0].std()
         background[i] = bkg
 
-    if np.abs(subtract_error) > 0:
-        n_data_array, n_error_array, background, error_bkg = subtract_bkg(n_data_array, n_error_array, mask, background, error_bkg)
-    
     if display:
         display_bkg(data, background, std_bkg, headers, rectangle=rectangle, savename=savename, plots_folder=plots_folder)
-    return n_data_array, n_error_array, headers, background, error_bkg
+    return n_data_array, n_error_array, headers, background
-
-def subtract_bkg(data, error, mask, background, error_bkg):
-    assert data.ndim == 3, "Input data must have more than 1 image."
-    
-    n_data_array, n_error_array = deepcopy(data), deepcopy(error)
-    
-    for i in range(data.shape[0]):
-        n_data_array[i][mask] = n_data_array[i][mask] - background[i]
-        n_data_array[i][np.logical_and(mask, n_data_array[i] <= 1e-3 * background[i])] = 1e-3 * background[i]
-        n_error_array[i] = np.sqrt(n_error_array[i]**2 + error_bkg[i]**2)
-    
-    return n_data_array, n_error_array, background, error_bkg

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 = [upsampled_region_size] * data.ndim
-            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 = [
+        axis_offsets = [0] * data.ndim
-            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
 

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/fits.py

@@ -16,7 +16,7 @@ from astropy.io import fits
 from astropy.wcs import WCS
 
 from .convex_hull import clean_ROI
-from .utils import wcs_PA
+from .utils import wcs_CD_to_PC, wcs_PA
 
 
 def get_obs_data(infiles, data_folder="", compute_flux=False):
@@ -46,40 +46,54 @@ def get_obs_data(infiles, data_folder="", compute_flux=False):
             data_array.append(f[0].data)
             wcs_array.append(WCS(header=f[0].header, fobj=f).celestial)
             f.flush()
+        # Save pixel area for flux density computation
+        if "PXFORMT" in headers[i].keys() and headers[i]["PXFORMT"] == "NORMAL":
+            headers[i]["PXAREA"] = 1.96e-4  # 14x14 milliarcsec squared pixel area in arcsec^2
+        elif "PXFORMT" in headers[i].keys() and headers[i]["PXFORMT"] == "ZOOM":
+            headers[i]["PXAREA"] = 4.06e-4  # 29x14 milliarcsec squared pixel area in arcsec^2
+        else:
+            headers[i]["PXAREA"] = 1.0  # unknown default to 1 arcsec^2
+        # Convert PHOTFLAM value from 1arcsec aperture to the pixel area
+        # headers[i]["PHOTFLAM"] *= np.pi / headers[i]["PXAREA"]
     data_array = np.array(data_array, dtype=np.double)
 
     # Prevent negative count value in imported data
     for i in range(len(data_array)):
         data_array[i][data_array[i] < 0.0] = 0.0
 
-    # force WCS to convention PCi_ja unitary, cdelt in deg
+    # Compute CDELT, ORIENTAT from header
     for wcs, header in zip(wcs_array, headers):
         new_wcs = wcs.deepcopy()
-        if new_wcs.wcs.has_cd() or (new_wcs.wcs.cdelt[:2] == np.array([1.0, 1.0])).all():
-            # Update WCS with relevant information
         if new_wcs.wcs.has_cd():
+            # Update WCS with relevant information
             del new_wcs.wcs.cd
             keys = list(new_wcs.to_header().keys()) + ["CD1_1", "CD1_2", "CD1_3", "CD2_1", "CD2_2", "CD2_3", "CD3_1", "CD3_2", "CD3_3"]
             for key in keys:
                 header.remove(key, ignore_missing=True)
             new_cdelt = np.linalg.eigvals(wcs.wcs.cd)
-                new_cdelt.sort()
             new_wcs.wcs.pc = wcs.wcs.cd.dot(np.diag(1.0 / new_cdelt))
             new_wcs.wcs.cdelt = new_cdelt
+            try:
+                header["ORIENTAT"] = float(header["ORIENTAT"])
+            except KeyError:
+                header["ORIENTAT"] = -wcs_CD_to_PC(new_wcs.wcs.cd)[1]
+        elif (new_wcs.wcs.cdelt[:2] != np.array([1.0, 1.0])).all():
+            try:
+                header["ORIENTAT"] = float(header["ORIENTAT"])
+            except KeyError:
+                header["ORIENTAT"] = -wcs_PA(new_wcs.wcs.pc, new_wcs.wcs.cdelt)
+        else:
+            print("Couldn't compute ORIENTAT from WCS")
         for key, val in new_wcs.to_header().items():
             header[key] = val
-        try:
-            _ = header["ORIENTAT"]
-        except KeyError:
-            header["ORIENTAT"] = wcs_PA(new_wcs.wcs.pc[1, 0], np.diag(new_wcs.wcs.pc).mean())
 
     # force WCS for POL60 to have same pixel size as POL0 and POL120
     is_pol60 = np.array([head["filtnam1"].lower() == "pol60" for head in headers], dtype=bool)
     cdelt = np.round(np.array([WCS(head).wcs.cdelt[:2] for head in headers]), 10)
-    if np.unique(cdelt[np.logical_not(is_pol60)], axis=0).size != 2:
+    if np.any(is_pol60) and np.unique(cdelt[np.logical_not(is_pol60)], axis=0).size != 2:
         print(np.unique(cdelt[np.logical_not(is_pol60)], axis=0))
         raise ValueError("Not all images have same pixel size")
-    else:
+    elif np.any(is_pol60):
         for i in np.arange(len(headers))[is_pol60]:
             headers[i]["cdelt1"], headers[i]["cdelt2"] = np.unique(cdelt[np.logical_not(is_pol60)], axis=0)[0]
 
@@ -92,7 +106,24 @@ def get_obs_data(infiles, data_folder="", compute_flux=False):
 
 
 def save_Stokes(
-    I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, header_stokes, data_mask, filename, data_folder="", return_hdul=False
+    I_stokes,
+    Q_stokes,
+    U_stokes,
+    Stokes_cov,
+    P,
+    debiased_P,
+    s_P,
+    s_P_P,
+    PA,
+    s_PA,
+    s_PA_P,
+    header_stokes,
+    data_mask,
+    filename,
+    data_folder="",
+    return_hdul=False,
+    flux_data=None,
+    flux_head=None,
 ):
     """
     Save computed polarimetry parameters to a single fits file,
@@ -141,12 +172,14 @@ def save_Stokes(
     header["TELESCOP"] = (header_stokes["TELESCOP"] if "TELESCOP" in list(header_stokes.keys()) else "HST", "telescope used to acquire data")
     header["INSTRUME"] = (header_stokes["INSTRUME"] if "INSTRUME" in list(header_stokes.keys()) else "FOC", "identifier for instrument used to acuire data")
     header["PHOTPLAM"] = (header_stokes["PHOTPLAM"], "Pivot Wavelength")
+    header["PHOTBW"] = (header_stokes["PHOTBW"], "RMS Bandwidth of the Filter and Detector")
     header["PHOTFLAM"] = (header_stokes["PHOTFLAM"], "Inverse Sensitivity in DN/sec/cm**2/Angst")
+    header["PXAREA"] = (header_stokes["PXAREA"], "Pixel area in arcsec**2")
     header["EXPTIME"] = (header_stokes["EXPTIME"], "Total exposure time in sec")
     header["PROPOSID"] = (header_stokes["PROPOSID"], "PEP proposal identifier for observation")
     header["TARGNAME"] = (header_stokes["TARGNAME"], "Target name")
     header["ORIENTAT"] = (header_stokes["ORIENTAT"], "Angle between North and the y-axis of the image")
-    header["FILENAME"] = (filename, "ORIGINAL FILENAME")
+    header["FILENAME"] = (filename, "Original filename")
     header["BKG_TYPE"] = (header_stokes["BKG_TYPE"], "Bkg estimation method used during reduction")
     header["BKG_SUB"] = (header_stokes["BKG_SUB"], "Amount of bkg subtracted from images")
     header["SMOOTH"] = (header_stokes["SMOOTH"] if "SMOOTH" in list(header_stokes.keys()) else "None", "Smoothing method used during reduction")
@@ -182,12 +215,30 @@ def save_Stokes(
     # Create HDUList object
     hdul = fits.HDUList([])
 
-    # Add I_stokes as PrimaryHDU
+    # Add Flux density as PrimaryHDU
-    header["datatype"] = ("I_stokes", "type of data stored in the HDU")
+    if flux_data is None:
+        header["DATATYPE"] = ("I_stokes", "type of data stored in the HDU")
         I_stokes[(1 - data_mask).astype(bool)] = 0.0
         primary_hdu = fits.PrimaryHDU(data=I_stokes, header=header)
         primary_hdu.name = "I_stokes"
         hdul.append(primary_hdu)
+    else:
+        flux_head["FILENAME"] = header["FILENAME"]
+        head = WCS(flux_head).deepcopy().to_header()
+        for key in [key for key in header.keys() if key not in ["SMOOTH", "SAMPLING"]]:
+            try:
+                head[key] = flux_head[key]
+            except KeyError:
+                head[key] = header[key]
+        header["DATATYPE"] = ("Flux_density", "type of data stored in the HDU")
+        primary_hdu = fits.PrimaryHDU(data=flux_data, header=head)
+        primary_hdu.name = "Flux_density"
+        hdul.append(primary_hdu)
+        header["DATATYPE"] = ("I_stokes", "type of data stored in the HDU")
+        I_stokes[(1 - data_mask).astype(bool)] = 0.0
+        image_hdu = fits.ImageHDU(data=I_stokes, header=header)
+        image_hdu.name = "I_stokes"
+        hdul.append(image_hdu)
 
     # Add Q, U, Stokes_cov, P, s_P, PA, s_PA to the HDUList
     for data, name in [
@@ -204,7 +255,7 @@ def save_Stokes(
         [data_mask, "Data_mask"],
     ]:
         hdu_header = header.copy()
-        hdu_header["datatype"] = name
+        hdu_header["DATATYPE"] = name
         if not name == "IQU_cov_matrix":
             data[(1 - data_mask).astype(bool)] = 0.0
         hdu = fits.ImageHDU(data=data, header=hdu_header)

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.
@@ -25,11 +25,11 @@ def divide_proposal(products):
     """
     for pid in np.unique(products["Proposal ID"]):
         obs = products[products["Proposal ID"] == pid].copy()
-        same_filt = np.unique(np.array(np.sum([obs["Filters"][:, 1:] == filt[1:] for filt in obs["Filters"]], axis=2) < 3, dtype=bool), axis=0)
+        same_filt = np.unique(np.array(np.sum([obs["Filters"] == filt for filt in obs["Filters"]], axis=2) >= len(obs["Filters"][0]), dtype=bool), axis=0)
         if len(same_filt) > 1:
             for filt in same_filt:
                 products["Proposal ID"][np.any([products["Dataset"] == dataset for dataset in obs["Dataset"][filt]], axis=0)] = "_".join(
-                    [obs["Proposal ID"][filt][0], "_".join([fi for fi in obs["Filters"][filt][0][1:] if fi[:-1] != "CLEAR"])]
+                    [obs["Proposal ID"][filt][0], "_".join([fi for fi in obs["Filters"][filt][0] if fi[:-1] != "CLEAR"])]
                 )
     for pid in np.unique(products["Proposal ID"]):
         obs = products[products["Proposal ID"] == pid].copy()
@@ -44,37 +44,68 @@ def divide_proposal(products):
     return products
 
 
-def get_product_list(target=None, proposal_id=None):
+def get_product_list(target=None, proposal_id=None, instrument="foc"):
     """
     Retrieve products list for a given target from the MAST archive
     """
     mission = MastMissions(mission="hst")
     radius = "3"
     select_cols = [
-        "sci_data_set_name",
-        "sci_spec_1234",
-        "sci_actual_duration",
-        "sci_start_time",
-        "sci_stop_time",
-        "sci_central_wavelength",
-        "sci_instrume",
-        "sci_aper_1234",
-        "sci_targname",
         "sci_pep_id",
         "sci_pi_last_name",
+        "sci_targname",
+        "sci_aper_1234",
+        "sci_spec_1234",
+        "sci_central_wavelength",
+        "sci_actual_duration",
+        "sci_instrume",
+        "sci_operating_mode",
+        "sci_data_set_name",
+        "sci_start_time",
+        "sci_stop_time",
+        "sci_refnum",
     ]
 
-    cols = ["Dataset", "Filters", "Exptime", "Start", "Stop", "Central wavelength", "Instrument", "Size", "Target name", "Proposal ID", "PI last name"]
+    cols = [
+        "Proposal ID",
+        "PI last name",
+        "Target name",
+        "Aperture",
+        "Filters",
+        "Central wavelength",
+        "Exptime",
+        "Instrument",
+        "Operating Mode",
+        "Dataset",
+        "Start",
+        "Stop",
+        "References",
+    ]
 
     if target is None:
         target = input("Target name:\n>")
 
     # Use query_object method to resolve the object name into coordinates
-    results = mission.query_object(target, radius=radius, select_cols=select_cols, sci_spec_1234="POL*", sci_obs_type="image", sci_aec="S", sci_instrume="foc")
+    if instrument == "foc":
+        results = mission.query_object(
+            target, radius=radius, select_cols=select_cols, sci_spec_1234="POL*", sci_obs_type="image", sci_aec="S", sci_instrume="foc"
+        )
+        dataproduct_type = "image"
+        description = "DADS C0F file - Calibrated exposure WFPC/WFPC2/FOC/FOS/GHRS/HSP"
+    elif instrument == "fos":
+        results = mission.query_object(
+            target, radius=radius, select_cols=select_cols, sci_operating_mode="SPECTROPOLARIMETRY", sci_obs_type="spectrum", sci_aec="S", sci_instrume="fos"
+        )
+        dataproduct_type = "spectrum"
+        description = ["DADS C0F file - Calibrated exposure WFPC/WFPC2/FOC/FOS/GHRS/HSP", "DADS C3F file - Calibrated exposure GHRS/FOS/HSP"]
 
     for c, n_c in zip(select_cols, cols):
         results.rename_column(c, n_c)
     results["Proposal ID"] = Column(results["Proposal ID"], dtype="U35")
+    if instrument == "foc":
+        results["POLFilters"] = Column(np.array([filt.split(";")[0] for filt in results["Filters"]], dtype=str))
+        results["Filters"] = Column(np.array([filt.split(";")[1:] for filt in results["Filters"]], dtype=str))
+    else:
         results["Filters"] = Column(np.array([filt.split(";") for filt in results["Filters"]], dtype=str))
     results["Start"] = Column(Time(results["Start"]))
     results["Stop"] = Column(Time(results["Stop"]))
@@ -89,20 +120,34 @@ def get_product_list(target=None, proposal_id=None):
             to_remove.append(i)
     obs.remove_rows(to_remove)
     # Remove observations for which a polarization filter is missing
+    if instrument == "foc":
         polfilt = {"POL0": 0, "POL60": 1, "POL120": 2}
         for pid in np.unique(obs["Proposal ID"]):
             used_pol = np.zeros(3)
             for dataset in obs[obs["Proposal ID"] == pid]:
-            used_pol[polfilt[dataset["Filters"][0]]] += 1
+                used_pol[polfilt[dataset["POLFilters"]]] += 1
             if np.any(used_pol < 1):
                 obs.remove_rows(np.arange(len(obs))[obs["Proposal ID"] == pid])
+    # Remove observations for which a spectropolarization has not been reduced
+    if instrument == "fos":
+        for pid in np.unique(obs["Proposal ID"]):
+            observations = Observations.query_criteria(proposal_id=pid.split("_")[0])
+            c3prod = Observations.filter_products(
+                Observations.get_product_list(observations),
+                productType=["SCIENCE"],
+                dataproduct_type=dataproduct_type,
+                calib_level=[2],
+                description=description[1],
+            )
+            if len(c3prod) < 1:
+                obs.remove_rows(np.arange(len(obs))[obs["Proposal ID"] == pid])
 
     tab = unique(obs, ["Target name", "Proposal ID"])
     obs["Obs"] = [np.argmax(np.logical_and(tab["Proposal ID"] == data["Proposal ID"], tab["Target name"] == data["Target name"])) + 1 for data in obs]
     try:
-        n_obs = unique(obs[["Obs", "Filters", "Start", "Central wavelength", "Instrument", "Size", "Target name", "Proposal ID", "PI last name"]], "Obs")
+        n_obs = unique(obs[["Obs", "Filters", "Start", "Central wavelength", "Instrument", "Aperture", "Target name", "Proposal ID", "PI last name"]], "Obs")
     except IndexError:
-        raise ValueError("There is no observation with POL0, POL60 and POL120 for {0:s} in HST/FOC Legacy Archive".format(target))
+        raise ValueError("There is no observation with polarimetry for {0:s} in HST/{1:s} Legacy Archive".format(target, instrument.upper()))
 
     b = np.zeros(len(results), dtype=bool)
     if proposal_id is not None and str(proposal_id) in obs["Proposal ID"]:
@@ -128,31 +173,25 @@ def get_product_list(target=None, proposal_id=None):
 
     observations = Observations.query_criteria(obs_id=list(results["Dataset"][b]))
     products = Observations.filter_products(
-        Observations.get_product_list(observations),
+        Observations.get_product_list(observations), productType=["SCIENCE"], dataproduct_type=dataproduct_type, calib_level=[2], description=description
-        productType=["SCIENCE"],
-        dataproduct_type=["image"],
-        calib_level=[2],
-        description="DADS C0F file - Calibrated exposure WFPC/WFPC2/FOC/FOS/GHRS/HSP",
     )
+
     products["proposal_id"] = Column(products["proposal_id"], dtype="U35")
-    products["target_name"] = Column(observations["target_name"])
 
     for prod in products:
         prod["proposal_id"] = results["Proposal ID"][results["Dataset"] == prod["productFilename"][: len(results["Dataset"][0])].upper()][0]
 
-    for prod in products:
+    tab = unique(products, "proposal_id")
-        prod["target_name"] = observations["target_name"][observations["obsid"] == prod["obsID"]][0]
-    tab = unique(products, ["target_name", "proposal_id"])
 
-    products["Obs"] = [np.argmax(np.logical_and(tab["proposal_id"] == data["proposal_id"], tab["target_name"] == data["target_name"])) + 1 for data in products]
+    products["Obs"] = [np.argmax(tab["proposal_id"] == data["proposal_id"]) + 1 for data in products]
     return target, products
 
 
-def retrieve_products(target=None, proposal_id=None, output_dir="./data"):
+def retrieve_products(target=None, proposal_id=None, instrument="foc", output_dir="./data"):
     """
     Given a target name and a proposal_id, create the local directories and retrieve the fits files from the MAST Archive
     """
-    target, products = get_product_list(target=target, proposal_id=proposal_id)
+    target, products = get_product_list(target=target, proposal_id=proposal_id, instrument=instrument)
     prodpaths = []
     # data_dir = path_join(output_dir, target)
     out = ""
@@ -183,10 +222,11 @@ if __name__ == "__main__":
     parser = argparse.ArgumentParser(description="Query MAST for target products")
     parser.add_argument("-t", "--target", metavar="targetname", required=False, help="the name of the target", type=str, default=None)
     parser.add_argument("-p", "--proposal_id", metavar="proposal_id", required=False, help="the proposal id of the data products", type=int, default=None)
+    parser.add_argument("-i", "--instrum", metavar="instrum", required=False, help="the instrument used for observation", type=str, default="foc")
     parser.add_argument(
         "-o", "--output_dir", metavar="directory_path", required=False, help="output directory path for the data products", type=str, default="./data"
     )
     args = parser.parse_args()
     print(args.target)
-    prodpaths = retrieve_products(target=args.target, proposal_id=args.proposal_id)
+    prodpaths = retrieve_products(target=args.target, proposal_id=args.proposal_id, instrument=args.instrum.lower(), output_dir=args.output_dir)
     print(prodpaths)

package/lib/reduction.py

@@ -52,7 +52,6 @@ from scipy.ndimage import rotate as sc_rotate
 from scipy.ndimage import shift as sc_shift
 from scipy.signal import fftconvolve
 
-from .background import bkg_fit, bkg_hist, bkg_mini
 from .convex_hull import image_hull
 from .cross_correlation import phase_cross_correlation
 from .deconvolve import deconvolve_im, gaussian2d, gaussian_psf, zeropad
@@ -192,7 +191,7 @@ 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')
+    >>> n = bin_ndarray(m, new_shape=(5, 5), operation="sum")
     >>> print(n)
 
     [[ 22  30  38  46  54]
@@ -279,9 +278,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 = [null_val] * error_array.shape[0]
-            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
@@ -310,6 +307,8 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu
         # Update CRPIX value in the associated header
         curr_wcs = WCS(crop_headers[i]).celestial.deepcopy()
         curr_wcs.wcs.crpix[:2] = curr_wcs.wcs.crpix[:2] - np.array([v_array[2], v_array[0]])
+        curr_wcs.array_shape = crop_array[i].shape
+        curr_wcs.wcs.set()
         crop_headers[i].update(curr_wcs.to_header())
         crop_headers[i]["naxis1"], crop_headers[i]["naxis2"] = crop_array[i].shape
 
@@ -348,10 +347,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=[rectangle] * len(headers),
-                    rectangle,
-                ]
-                * len(headers),
                 savename=savename + "_crop_region",
                 plots_folder=plots_folder,
             )
@@ -416,12 +412,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))
 
@@ -446,9 +442,9 @@ def get_error(
     return_background=False,
 ):
     """
-    Look for sub-image of shape sub_shape that have the smallest integrated
+    Estimate background intensity level from either fitting the intensity histogram
-    flux (no source assumption) and define the background on the image by the
+    or by looking for the sub-image of smallest integrated intensity (no source assumption)
-    standard deviation on this sub-image.
+    and define the background on the image by the standard deviation on this sub-image.
     ----------
     Inputs:
     data_array : numpy.ndarray
@@ -468,7 +464,7 @@ def get_error(
         If 'auto', look for optimal binning and fit intensity histogram with au gaussian.
         If str or None, statistic rule to be used for the number of bins in counts/s.
         If int, number of bins for the counts/s histogram.
-        If tuple, shape of the sub-image to look for. Must be odd.
+        If tuple, shape of the sub-image of lowest intensity to look for.
         Defaults to None.
     subtract_error : float or bool, optional
         If float, factor to which the estimated background should be multiplied
@@ -528,25 +524,29 @@ def get_error(
     # estimated to less than 3%
     err_flat = data * 0.03
 
+    from .background import bkg_fit, bkg_hist, bkg_mini
 
-    if (sub_type is None):
+    if sub_type is None:
-        n_data_array, c_error_bkg, headers, background, error_bkg = bkg_hist(
+        n_data_array, c_error_bkg, headers, background = bkg_hist(
-            data, error, mask, headers, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder)
+            data, error, mask, headers, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder
+        )
         sub_type, subtract_error = "histogram ", str(int(subtract_error > 0.0))
     elif isinstance(sub_type, str):
-        if sub_type.lower() in ['auto']:
+        if sub_type.lower() in ["fit"]:
-            n_data_array, c_error_bkg, headers, background, error_bkg = bkg_fit(
+            n_data_array, c_error_bkg, headers, background = bkg_fit(
-                data, error, mask, headers, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder)
+                data, error, mask, headers, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder
+            )
             sub_type, subtract_error = "histogram fit ", "mean+%.1fsigma" % subtract_error if subtract_error != 0.0 else 0.0
         else:
-            n_data_array, c_error_bkg, headers, background, error_bkg = bkg_hist(
+            n_data_array, c_error_bkg, headers, background = bkg_hist(
-                data, error, mask, headers, sub_type=sub_type, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder)
+                data, error, mask, headers, n_bins=sub_type, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder
+            )
             sub_type, subtract_error = "histogram ", "mean+%.1fsigma" % subtract_error if subtract_error != 0.0 else 0.0
     elif isinstance(sub_type, tuple):
-        n_data_array, c_error_bkg, headers, background, error_bkg = bkg_mini(
+        n_data_array, c_error_bkg, headers, background = bkg_mini(
-            data, error, mask, headers, sub_shape=sub_type, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder)
+            data, error, mask, headers, sub_shape=sub_type, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder
+        )
         sub_type, subtract_error = "minimal flux ", "mean+%.1fsigma" % subtract_error if subtract_error != 0.0 else 0.0
-
     else:
         print("Warning: Invalid subtype.")
 
@@ -558,7 +558,7 @@ def get_error(
     n_error_array = np.sqrt(err_wav**2 + err_psf**2 + err_flat**2 + c_error_bkg**2)
 
     if return_background:
-        return n_data_array, n_error_array, headers, background, error_bkg # return background error as well
+        return n_data_array, n_error_array, headers, background
     else:
         return n_data_array, n_error_array, headers
 
@@ -627,12 +627,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(
+            Dxy_arr[i] = np.array([pxsize] * 2)
-                [
-                    pxsize,
-                ]
-                * 2
-            )
             scale = "px"
         elif scale.lower() in ["arcsec", "arcseconds"]:
             Dxy_arr[i] = np.array(pxsize / np.abs(w.wcs.cdelt) / 3600.0)
@@ -642,7 +637,8 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio
             pxsize, scale = "", "full"
         else:
             raise ValueError("'{0:s}' invalid scale for binning.".format(scale))
-    new_shape = np.ceil(min(image.shape / Dxy_arr, key=lambda x: x[0] + x[1])).astype(int)
+    new_shape_float = min(image.shape / Dxy_arr, key=lambda x: x[0] + x[1])
+    new_shape = np.ceil(new_shape_float).astype(int)
 
     for i, (image, error, header) in enumerate(list(zip(data_array, error_array, headers))):
         # Get current pixel size
@@ -671,9 +667,12 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio
         # Update header
         nw = w.deepcopy()
         nw.wcs.cdelt *= Dxy
+        # nw.wcs.crpix += np.abs(new_shape_float - new_shape) * np.array(new_shape) / Dxy
         nw.wcs.crpix /= Dxy
         nw.array_shape = new_shape
+        nw.wcs.set()
         new_header["NAXIS1"], new_header["NAXIS2"] = nw.array_shape
+        new_header["PXAREA"] *= Dxy[0] * Dxy[1]
         for key, val in nw.to_header().items():
             new_header.set(key, val)
         new_header["SAMPLING"] = (str(pxsize) + scale, "Resampling performed during reduction")
@@ -691,7 +690,7 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio
 
 
 def align_data(
-    data_array, headers, error_array=None, data_mask=None, background=None, upsample_factor=1.0, ref_data=None, ref_center=None, return_shifts=False, optimal_binning=False
+    data_array, headers, error_array=None, data_mask=None, background=None, upsample_factor=1.0, ref_data=None, ref_center=None, return_shifts=False
 ):
     """
     Align images in data_array using cross correlation, and rescale them to
@@ -770,7 +769,6 @@ def align_data(
     full_headers.append(headers[0])
     err_array = np.concatenate((error_array, [np.zeros(ref_data.shape)]), axis=0)
 
-    if not optimal_binning:
     if data_mask is None:
         full_array, err_array, full_headers = crop_array(full_array, full_headers, err_array, step=5, inside=False, null_val=0.0)
     else:
@@ -852,11 +850,12 @@ def align_data(
     new_crpix = np.array([wcs.wcs.crpix for wcs in headers_wcs]) + shifts[:, ::-1] + res_shift[::-1]
     for i in range(len(headers_wcs)):
         headers_wcs[i].wcs.crpix = new_crpix[0]
+        headers_wcs[i].array_shape = (res_shape, res_shape)
+        headers_wcs[i].wcs.set()
         headers[i].update(headers_wcs[i].to_header())
+        headers[i]["NAXIS1"], headers[i]["NAXIS2"] = res_shape, res_shape
 
     data_mask = rescaled_mask.all(axis=0)
-    
-    if not optimal_binning:
     data_array, error_array, data_mask, headers = crop_array(rescaled_image, headers, rescaled_error, data_mask, null_val=0.01 * background)
 
     if return_shifts:
@@ -938,12 +937,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(
+                    g_rc = np.array([np.exp(-0.5 * (dist_rc / stdev) ** 2) / (2.0 * np.pi * stdev**2)] * data_array.shape[0])
-                        [
-                            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 +1432,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 = [
+        all_header_stokes = [{}] * np.unique(rotate).size
-            {},
-        ] * np.unique(rotate).size
 
         for i, rot in enumerate(np.unique(rotate)):
             rot_mask = rotate == rot
@@ -1723,13 +1715,11 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_st
 
     new_wcs.wcs.pc = np.dot(mrot, new_wcs.wcs.pc)
     new_wcs.wcs.crpix = np.dot(mrot, new_wcs.wcs.crpix - old_center[::-1]) + new_center[::-1]
+    new_wcs.array_shape = shape
     new_wcs.wcs.set()
     for key, val in new_wcs.to_header().items():
         new_header_stokes.set(key, val)
-    if new_wcs.wcs.pc[0, 0] == 1.0:
+    new_header_stokes["NAXIS1"], new_header_stokes["NAXIS2"] = new_wcs.array_shape
-        new_header_stokes.set("PC1_1", 1.0)
-    if new_wcs.wcs.pc[1, 1] == 1.0:
-        new_header_stokes.set("PC2_2", 1.0)
     new_header_stokes["ORIENTAT"] += ang
 
     # Nan handling :
@@ -1802,8 +1792,6 @@ def rotate_data(data_array, error_array, data_mask, headers):
         Updated list of headers corresponding to the reduced images accounting
         for the new orientation angle.
     """
-    # Rotate I_stokes, Q_stokes, U_stokes using rotation matrix
-
     old_center = np.array(data_array[0].shape) / 2
     shape = np.fix(np.array(data_array[0].shape) * np.sqrt(2.5)).astype(int)
     new_center = np.array(shape) / 2
@@ -1832,9 +1820,11 @@ def rotate_data(data_array, error_array, data_mask, headers):
         new_wcs = WCS(header).celestial.deepcopy()
         new_wcs.wcs.pc[:2, :2] = np.dot(mrot, new_wcs.wcs.pc[:2, :2])
         new_wcs.wcs.crpix[:2] = np.dot(mrot, new_wcs.wcs.crpix[:2] - old_center[::-1]) + new_center[::-1]
+        new_wcs.array_shape = shape
         new_wcs.wcs.set()
         for key, val in new_wcs.to_header().items():
             new_header[key] = val
+        new_header["NAXIS1"], new_header["NAXIS2"] = new_wcs.array_shape
         new_header["ORIENTAT"] = np.arccos(new_wcs.celestial.wcs.pc[0, 0]) * 180.0 / np.pi
         new_header["ROTATE"] = ang
         new_headers.append(new_header)

package/lib/utils.py

@@ -4,6 +4,14 @@ import numpy as np
 def rot2D(ang):
     """
     Return the 2D rotation matrix of given angle in degrees
+    ----------
+    Inputs:
+    ang : float
+        Angle in degrees
+    ----------
+    Returns:
+    rot_mat : numpy.ndarray
+        2D matrix of shape (2,2) to perform vector rotation at angle "ang".
     """
     alpha = np.pi * ang / 180
     return np.array([[np.cos(alpha), np.sin(alpha)], [-np.sin(alpha), np.cos(alpha)]])
@@ -13,6 +21,14 @@ def princ_angle(ang):
     """
     Return the principal angle in the 0° to 180° quadrant as PA is always
     defined at p/m 180°.
+    ----------
+    Inputs:
+    ang : float, numpy.ndarray
+        Angle in degrees. Can be an array of angles.
+    ----------
+    Returns:
+    princ_ang : float, numpy.ndarray
+        Principal angle in the 0°-180° quadrant in the same shape as input.
     """
     if not isinstance(ang, np.ndarray):
         A = np.array([ang])
@@ -28,9 +44,100 @@ def princ_angle(ang):
         return A[0]
 
 
+def PCconf(QN, UN, QN_ERR, UN_ERR):
+    """
+    Compute the confidence level for 2 parameters polarisation degree and
+    polarisation angle from the PCUBE analysis.
+    ----------
+    Inputs:
+    QN : float, numpy.ndarray
+        Normalized Q Stokes flux.
+    UN : float, numpy.ndarray
+        Normalized U Stokes flux.
+    QN_ERR : float, numpy.ndarray
+        Normalized error on Q Stokes flux.
+    UN_ERR : float, numpy.ndarray
+        Normalized error on U Stokes flux.
+    ----------
+    Returns:
+    conf : numpy.ndarray
+        2D matrix of same shape as input containing the confidence on the polarization
+        computation between 0 and 1 for 2 parameters of interest (Q and U Stokes fluxes).
+    """
+    mask = np.logical_and(QN_ERR > 0.0, UN_ERR > 0.0)
+    conf = np.full(QN.shape, -1.0)
+    chi2 = QN**2 / QN_ERR**2 + UN**2 / UN_ERR**2
+    conf[mask] = 1.0 - np.exp(-0.5 * chi2[mask])
+    return conf
+
+
+def CenterConf(mask, PA, sPA):
+    """
+    Compute the confidence map for the position of the center of emission.
+    ----------
+    Inputs:
+    mask : bool, numpy.ndarray
+        Mask of the polarization vectors from which the center of emission should be drawn.
+    PA : float, numpy.ndarray
+        2D matrix containing the computed polarization angle.
+    sPA : float, numpy.ndarray
+        2D matrix containing the total uncertainties on the polarization angle.
+    ----------
+    Returns:
+    conf : numpy.ndarray
+        2D matrix of same shape as input containing the confidence on the polarization
+        computation between 0 and 1 for 2 parameters of interest (Q and U Stokes fluxes).
+    """
+    chi2 = np.full(PA.shape, np.nan, dtype=np.float64)
+    conf = np.full(PA.shape, -1.0, dtype=np.float64)
+    yy, xx = np.indices(PA.shape)
+    Nobs = np.sum(mask)
+
+    def ideal(c):
+        itheta = np.full(PA.shape, np.nan)
+        itheta[(xx + 0.5) != c[0]] = np.degrees(np.arctan((yy[(xx + 0.5) != c[0]] + 0.5 - c[1]) / (xx[(xx + 0.5) != c[0]] + 0.5 - c[0])))
+        itheta[(xx + 0.5) == c[0]] = PA[(xx + 0.5) == c[0]]
+        return princ_angle(itheta)
+
+    def chisq(c):
+        return np.sum((princ_angle(PA[mask]) - ideal((c[0], c[1]))[mask]) ** 2 / sPA[mask] ** 2) / (Nobs - 2)
+
+    for x, y in zip(xx[np.isfinite(PA)], yy[np.isfinite(PA)]):
+        chi2[y, x] = chisq((x, y))
+
+    from scipy.optimize import minimize
+
+    conf[np.isfinite(PA)] = 0.5 * np.exp(-0.5 * chi2[np.isfinite(PA)])
+    c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1]
+    result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr")
+    if result.success:
+        print("Center of emission found: reduced chi_squared {0:.2f}/{1:d}={2:.2f}".format(chisq(result.x) * (Nobs - 2), Nobs - 2, chisq(result.x)))
+    else:
+        print("Center of emission not found", result)
+    return conf, result.x
+
+
 def sci_not(v, err, rnd=1, out=str):
     """
-    Return the scientifque error notation as a string.
+    Return the scientific error notation as a string.
+    ----------
+    Inputs:
+    v : float
+        Value to be transformed into scientific notation.
+    err : float
+        Error on the value to be transformed into scientific notation.
+    rnd : int
+        Number of significant numbers for the scientific notation.
+        Default to 1.
+    out : str or other
+        Format in which the notation should be returned. "str" means the notation
+        is returned as a single string, "other" means it is returned as a list of "str".
+        Default to str.
+    ----------
+    Returns:
+    conf : numpy.ndarray
+        2D matrix of same shape as input containing the confidence on the polarization
+        computation between 0 and 1 for 2 parameters of interest (Q and U Stokes fluxes).
     """
     power = -int(("%E" % v)[-3:]) + 1
     output = [r"({0}".format(round(v * 10**power, rnd)), round(v * 10**power, rnd)]
@@ -46,17 +153,47 @@ def sci_not(v, err, rnd=1, out=str):
     else:
         return *output[1:], -power
 
-def wcs_PA(PC21, PC22):
+
+def wcs_CD_to_PC(CD):
     """
     Return the position angle in degrees to the North direction of a wcs
     from the values of coefficient of its transformation matrix.
+    ----------
+    Inputs:
+    CD : np.ndarray
+        Value of the WCS matrix CD
+    ----------
+    Returns:
+    cdelt : (float, float)
+        Scaling factor in arcsec between pixel in X and Y directions.
+    orient : float
+        Angle in degrees between the North direction and the Up direction of the WCS.
     """
-    if (abs(PC21) > abs(PC22)) and (PC21 >= 0):
+    det = CD[0, 0] * CD[1, 1] - CD[0, 1] * CD[1, 0]
-        orient = -np.arccos(PC22) * 180.0 / np.pi
+    sgn = -1.0 if det < 0 else 1.0
-    elif (abs(PC21) > abs(PC22)) and (PC21 < 0):
+    cdelt = np.array([sgn, 1.0]) * np.sqrt(np.sum(CD**2, axis=1))
-        orient = np.arccos(PC22) * 180.0 / np.pi
+    rot = np.arctan2(-CD[1, 0], sgn * CD[0, 0])
-    elif (abs(PC21) < abs(PC22)) and (PC22 >= 0):
+    rot2 = np.arctan2(sgn * CD[0, 1], CD[1, 1])
-        orient = np.arccos(PC22) * 180.0 / np.pi
+    orient = 0.5 * (rot + rot2) * 180.0 / np.pi
-    elif (abs(PC21) < abs(PC22)) and (PC22 < 0):
+    return cdelt, orient
-        orient = -np.arccos(PC22) * 180.0 / np.pi
+
+
+def wcs_PA(PC, cdelt):
+    """
+    Return the position angle in degrees to the North direction of a wcs
+    from the values of coefficient of its transformation matrix.
+    ----------
+    Inputs:
+    PC : np.ndarray
+        Value of the WCS matrix PC
+    cdelt : (float, float)
+        Scaling factor in arcsec between pixel in X and Y directions.
+    ----------
+    Returns:
+    orient : float
+        Angle in degrees between the North direction and the Up direction of the WCS.
+    """
+    rot = np.pi / 2.0 - np.arctan2(-cdelt[1] * PC[1, 0], abs(cdelt[0]) * PC[0, 0])
+    rot2 = np.pi / 2.0 - np.arctan2(abs(cdelt[0]) * PC[0, 1], cdelt[1] * PC[1, 1])
+    orient = 0.5 * (rot + rot2) * 180.0 / np.pi
     return orient

package/overplot_IC5063.py

@@ -1,50 +0,0 @@
-#!/usr/bin/python3
-import numpy as np
-from astropy.io import fits
-from lib.plots import overplot_pol, overplot_radio
-from matplotlib.colors import LogNorm
-
-Stokes_UV = fits.open("./data/IC5063/5918/IC5063_FOC_b0.10arcsec_c0.20arcsec.fits")
-Stokes_18GHz = fits.open("./data/IC5063/radio/IC5063_18GHz.fits")
-Stokes_24GHz = fits.open("./data/IC5063/radio/IC5063_24GHz.fits")
-Stokes_103GHz = fits.open("./data/IC5063/radio/IC5063_103GHz.fits")
-Stokes_229GHz = fits.open("./data/IC5063/radio/IC5063_229GHz.fits")
-Stokes_357GHz = fits.open("./data/IC5063/radio/IC5063_357GHz.fits")
-# Stokes_S2 = fits.open("./data/IC5063/POLARIZATION_COMPARISON/S2_rot_crop.fits")
-Stokes_IR = fits.open("./data/IC5063/IR/u2e65g01t_c0f_rot.fits")
-
-# levelsMorganti = np.array([1.,2.,3.,8.,16.,32.,64.,128.])
-levelsMorganti = np.logspace(-0.1249, 1.97, 7) / 100.0
-
-levels18GHz = levelsMorganti * Stokes_18GHz[0].data.max()
-A = overplot_radio(Stokes_UV, Stokes_18GHz)
-A.plot(levels=levels18GHz, SNRp_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/18GHz_overplot.pdf", vec_scale=None)
-
-levels24GHz = levelsMorganti * Stokes_24GHz[0].data.max()
-B = overplot_radio(Stokes_UV, Stokes_24GHz)
-B.plot(levels=levels24GHz, SNRp_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/24GHz_overplot.pdf", vec_scale=None)
-
-levels103GHz = levelsMorganti * Stokes_103GHz[0].data.max()
-C = overplot_radio(Stokes_UV, Stokes_103GHz)
-C.plot(levels=levels103GHz, SNRp_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/103GHz_overplot.pdf", vec_scale=None)
-
-levels229GHz = levelsMorganti * Stokes_229GHz[0].data.max()
-D = overplot_radio(Stokes_UV, Stokes_229GHz)
-D.plot(levels=levels229GHz, SNRp_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/229GHz_overplot.pdf", vec_scale=None)
-
-levels357GHz = levelsMorganti * Stokes_357GHz[0].data.max()
-E = overplot_radio(Stokes_UV, Stokes_357GHz)
-E.plot(levels=levels357GHz, SNRp_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/357GHz_overplot.pdf", vec_scale=None)
-
-# F = overplot_pol(Stokes_UV, Stokes_S2)
-# F.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='./plots/IC5063/S2_overplot.pdf', norm=LogNorm(vmin=5e-20,vmax=5e-18))
-
-G = overplot_pol(Stokes_UV, Stokes_IR, cmap="inferno")
-G.plot(
-    SNRp_cut=2.0,
-    SNRi_cut=10.0,
-    savename="./plots/IC5063/IR_overplot.pdf",
-    vec_scale=None,
-    norm=LogNorm(Stokes_IR[0].data.max() * Stokes_IR[0].header["photflam"] / 1e3, Stokes_IR[0].data.max() * Stokes_IR[0].header["photflam"]),
-    cmap="inferno_r",
-)

package/overplot_MRK463E.py

@@ -1,20 +0,0 @@
-#!/usr/bin/python3
-import numpy as np
-from astropy.io import fits
-from lib.plots import overplot_chandra, overplot_pol
-from matplotlib.colors import LogNorm
-
-Stokes_UV = fits.open("./data/MRK463E/5960/MRK463E_FOC_b0.05arcsec_c0.10arcsec.fits")
-Stokes_IR = fits.open("./data/MRK463E/WFPC2/IR_rot_crop.fits")
-Stokes_Xr = fits.open("./data/MRK463E/Chandra/X_ray_crop.fits")
-
-levels = np.geomspace(1.0, 99.0, 7)
-
-A = overplot_chandra(Stokes_UV, Stokes_Xr, norm=LogNorm())
-A.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=5, zoom=1, savename="./plots/MRK463E/Chandra_overplot.pdf")
-A.write_to(path1="./data/MRK463E/FOC_data_Chandra.fits", path2="./data/MRK463E/Chandra_data.fits", suffix="aligned")
-
-levels = np.array([0.8, 2, 5, 10, 20, 50]) / 100.0 * Stokes_UV[0].header["photflam"]
-B = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm())
-B.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=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")

package/src/Combine.py

@@ -1,6 +1,9 @@
-#!/usr/bin/python
+#!/usr/bin/env python3
 # -*- coding:utf-8 -*-
-# Project libraries
+from pathlib import Path
+from sys import path as syspath
+
+syspath.append(str(Path(__file__).parent.parent))
 
 import numpy as np
 
@@ -170,7 +173,7 @@ def main(infiles, target=None, output_dir="./data/"):
     # Reduction parameters
     kwargs = {}
     #  Polarization map output
-    kwargs["SNRp_cut"] = 3.0
+    kwargs["P_cut"] = 0.99
     kwargs["SNRi_cut"] = 1.0
     kwargs["flux_lim"] = 1e-19, 3e-17
     kwargs["scale_vec"] = 5
@@ -183,9 +186,7 @@ def main(infiles, target=None, output_dir="./data/"):
 
         new_infiles = []
         for i, group in enumerate(grouped_infiles):
-            new_infiles.append(
+            new_infiles.append(FOC_reduction(target=target + "-" + str(i + 1), infiles=["/".join([data_folder, file]) for file in group], interactive=True)[0])
-                FOC_reduction(target=target + "-" + str(i + 1), infiles=["/".join([data_folder, file]) for file in group], interactive=True)[0]
-            )
 
         infiles = new_infiles
 

package/src/analysis.py

@@ -1,40 +0,0 @@
-#!/usr/bin/python
-from getopt import error as get_error
-from getopt import getopt
-from sys import argv
-
-arglist = argv[1:]
-options = "hf:p:i:l:"
-long_options = ["help", "fits=", "snrp=", "snri=", "lim="]
-
-fits_path = None
-SNRp_cut, SNRi_cut = 3, 3
-flux_lim = None
-out_txt = None
-
-try:
-    arg, val = getopt(arglist, options, long_options)
-
-    for curr_arg, curr_val in arg:
-        if curr_arg in ("-h", "--help"):
-            print("python3 analysis.py -f <path_to_reduced_fits> -p <SNRp_cut> -i <SNRi_cut> -l <flux_lim>")
-        elif curr_arg in ("-f", "--fits"):
-            fits_path = str(curr_val)
-        elif curr_arg in ("-p", "--snrp"):
-            SNRp_cut = int(curr_val)
-        elif curr_arg in ("-i", "--snri"):
-            SNRi_cut = int(curr_val)
-        elif curr_arg in ("-l", "--lim"):
-            flux_lim = list("".join(curr_val).split(","))
-except get_error as err:
-    print(str(err))
-
-if fits_path is not None:
-    from astropy.io import fits
-    from lib.plots import pol_map
-
-    Stokes_UV = fits.open(fits_path)
-    p = pol_map(Stokes_UV, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim)
-
-else:
-    print("python3 analysis.py -f <path_to_reduced_fits> -p <SNRp_cut> -i <SNRi_cut> -l <flux_lim>")

package/src/comparison_Kishimoto.py

@@ -1,214 +0,0 @@
-#!/usr/bin/python
-from src.lib.background import gauss, bin_centers
-from src.lib.deconvolve import zeropad
-from src.lib.reduction import align_data
-from src.lib.plots import princ_angle
-from matplotlib.colors import LogNorm
-from os.path import join as path_join
-from astropy.io import fits
-from astropy.wcs import WCS
-from scipy.ndimage import shift
-from scipy.optimize import curve_fit
-import numpy as np
-import matplotlib.pyplot as plt
-
-root_dir = path_join('/home/t.barnouin/Documents/Thesis/HST')
-root_dir_K = path_join(root_dir, 'Kishimoto', 'output')
-root_dir_S = path_join(root_dir, 'FOC_Reduction', 'output')
-root_dir_data_S = path_join(root_dir, 'FOC_Reduction', 'data', 'NGC1068', '5144')
-root_dir_plot_S = path_join(root_dir, 'FOC_Reduction', 'plots', 'NGC1068', '5144', 'notaligned')
-filename_S = "NGC1068_FOC_b10.00pixel_not_aligned.fits"
-plt.rcParams.update({'font.size': 15})
-
-SNRi_cut = 30.
-SNRp_cut = 3.
-
-data_K = {}
-data_S = {}
-for d, i in zip(['I', 'Q', 'U', 'P', 'PA', 'sI', 'sQ', 'sU', 'sP', 'sPA'], [0, 1, 2, 5, 8, (3, 0, 0), (3, 1, 1), (3, 2, 2), 6, 9]):
-    data_K[d] = np.loadtxt(path_join(root_dir_K, d+'.txt'))
-    with fits.open(path_join(root_dir_data_S, filename_S)) as f:
-        if not type(i) is int:
-            data_S[d] = np.sqrt(f[i[0]].data[i[1], i[2]])
-        else:
-            data_S[d] = f[i].data
-    if i == 0:
-        header = f[i].header
-wcs = WCS(header)
-convert_flux = header['photflam']
-
-bkg_S = np.median(data_S['I'])/3
-bkg_K = np.median(data_K['I'])/3
-
-# zeropad data to get same size of array
-shape = data_S['I'].shape
-for d in data_K:
-    data_K[d] = zeropad(data_K[d], shape)
-
-# shift array to get same information in same pixel
-data_arr, error_ar, heads, data_msk, shifts, shifts_err = align_data(np.array([data_S['I'], data_K['I']]), [header, header], error_array=np.array(
-    [data_S['sI'], data_K['sI']]), background=np.array([bkg_S, bkg_K]), upsample_factor=10., ref_center='center', return_shifts=True)
-for d in data_K:
-    data_K[d] = shift(data_K[d], shifts[1], order=1, cval=0.)
-
-# compute pol components from shifted array
-for d in [data_S, data_K]:
-    for i in d:
-        d[i][np.isnan(d[i])] = 0.
-    d['P'] = np.where(np.logical_and(np.isfinite(d['I']), d['I'] > 0.), np.sqrt(d['Q']**2+d['U']**2)/d['I'], 0.)
-    d['sP'] = np.where(np.logical_and(np.isfinite(d['I']), d['I'] > 0.), np.sqrt((d['Q']**2*d['sQ']**2+d['U']**2*d['sU']**2) /
-                       (d['Q']**2+d['U']**2)+((d['Q']/d['I'])**2+(d['U']/d['I'])**2)*d['sI']**2)/d['I'], 0.)
-    d['d_P'] = np.where(np.logical_and(np.isfinite(d['P']), np.isfinite(d['sP'])), np.sqrt(d['P']**2-d['sP']**2), 0.)
-    d['PA'] = 0.5*np.arctan2(d['U'], d['Q'])+np.pi
-    d['SNRp'] = np.zeros(d['d_P'].shape)
-    d['SNRp'][d['sP'] > 0.] = d['d_P'][d['sP'] > 0.]/d['sP'][d['sP'] > 0.]
-    d['SNRi'] = np.zeros(d['I'].shape)
-    d['SNRi'][d['sI'] > 0.] = d['I'][d['sI'] > 0.]/d['sI'][d['sI'] > 0.]
-    d['mask'] = np.logical_and(d['SNRi'] > SNRi_cut, d['SNRp'] > SNRp_cut)
-data_S['mask'], data_K['mask'] = np.logical_and(data_S['mask'], data_K['mask']), np.logical_and(data_S['mask'], data_K['mask'])
-
-
-#
-#  Compute histogram of measured polarization in cut
-#
-bins = int(data_S['mask'].sum()/5)
-bin_size = 1./bins
-mod_p = np.linspace(0., 1., 300)
-for d in [data_S, data_K]:
-    d['hist'], d['bin_edges'] = np.histogram(d['d_P'][d['mask']], bins=bins, range=(0., 1.))
-    d['binning'] = bin_centers(d['bin_edges'])
-    peak, bins_fwhm = d['binning'][np.argmax(d['hist'])], d['binning'][d['hist'] > d['hist'].max()/2.]
-    fwhm = bins_fwhm[1]-bins_fwhm[0]
-    p0 = [d['hist'].max(), peak, fwhm]
-    try:
-        popt, pcov = curve_fit(gauss, d['binning'], d['hist'], p0=p0)
-    except RuntimeError:
-        popt = p0
-    d['hist_chi2'] = np.sum((d['hist']-gauss(d['binning'], *popt))**2)/d['hist'].size
-    d['hist_popt'] = popt
-
-fig_p, ax_p = plt.subplots(num="Polarization degree histogram", figsize=(10, 6), constrained_layout=True)
-ax_p.errorbar(data_S['binning'], data_S['hist'], xerr=bin_size/2., fmt='b.', ecolor='b', label='P through this pipeline')
-ax_p.plot(mod_p, gauss(mod_p, *data_S['hist_popt']), 'b--', label='mean = {1:.2f}, stdev = {2:.2f}'.format(*data_S['hist_popt']))
-ax_p.errorbar(data_K['binning'], data_K['hist'], xerr=bin_size/2., fmt='r.', ecolor='r', label="P through Kishimoto's pipeline")
-ax_p.plot(mod_p, gauss(mod_p, *data_K['hist_popt']), 'r--', label='mean = {1:.2f}, stdev = {2:.2f}'.format(*data_K['hist_popt']))
-ax_p.set(xlabel="Polarization degree", ylabel="Counts", title="Histogram of polarization degree computed in the cut for both pipelines.")
-ax_p.legend()
-fig_p.savefig(path_join(root_dir_plot_S, "NGC1068_K_pol_deg.png"), bbox_inches="tight", dpi=300)
-
-#
-#  Compute angular difference between the maps in cut
-#
-dtheta = np.where(data_S['mask'], 0.5*np.arctan((np.sin(2*data_S['PA'])*np.cos(2*data_K['PA'])-np.cos(2*data_S['PA']) *
-                  np.cos(2*data_K['PA']))/(np.cos(2*data_S['PA'])*np.cos(2*data_K['PA'])+np.cos(2*data_S['PA'])*np.sin(2*data_K['PA']))), np.nan)
-fig_pa = plt.figure(num="Polarization degree alignement")
-ax_pa = fig_pa.add_subplot(111, projection=wcs)
-cbar_ax_pa = fig_pa.add_axes([0.88, 0.12, 0.01, 0.75])
-ax_pa.set_title(r"Degree of alignement $\zeta$ of the polarization angles from the 2 pipelines in the cut")
-im_pa = ax_pa.imshow(np.cos(2*dtheta), vmin=-1., vmax=1., origin='lower', cmap='bwr', label=r"$\zeta$ between this pipeline and Kishimoto's")
-cbar_pa = plt.colorbar(im_pa, cax=cbar_ax_pa, label=r"$\zeta = \cos\left( 2 \cdot \delta\theta_P \right)$")
-ax_pa.coords[0].set_axislabel('Right Ascension (J2000)')
-ax_pa.coords[1].set_axislabel('Declination (J2000)')
-fig_pa.savefig(path_join(root_dir_plot_S, "NGC1068_K_pol_ang.png"), bbox_inches="tight", dpi=300)
-
-#
-#  Compute power uncertainty difference between the maps in cut
-#
-eta = np.where(data_S['mask'], np.abs(data_K['d_P']-data_S['d_P'])/np.sqrt(data_S['sP']**2+data_K['sP']**2)/2., np.nan)
-fig_dif_p = plt.figure(num="Polarization power difference ratio")
-ax_dif_p = fig_dif_p.add_subplot(111, projection=wcs)
-cbar_ax_dif_p = fig_dif_p.add_axes([0.88, 0.12, 0.01, 0.75])
-ax_dif_p.set_title(r"Degree of difference $\eta$ of the polarization from the 2 pipelines in the cut")
-im_dif_p = ax_dif_p.imshow(eta, vmin=0., vmax=2., origin='lower', cmap='bwr_r', label=r"$\eta$ between this pipeline and Kishimoto's")
-cbar_dif_p = plt.colorbar(im_dif_p, cax=cbar_ax_dif_p, label=r"$\eta = \frac{2 \left|P^K-P^S\right|}{\sqrt{{\sigma^K_P}^2+{\sigma^S_P}^2}}$")
-ax_dif_p.coords[0].set_axislabel('Right Ascension (J2000)')
-ax_dif_p.coords[1].set_axislabel('Declination (J2000)')
-fig_dif_p.savefig(path_join(root_dir_plot_S, "NGC1068_K_pol_diff.png"), bbox_inches="tight", dpi=300)
-
-#
-#  Compute angle uncertainty difference between the maps in cut
-#
-eta = np.where(data_S['mask'], np.abs(data_K['PA']-data_S['PA'])/np.sqrt(data_S['sPA']**2+data_K['sPA']**2)/2., np.nan)
-fig_dif_pa = plt.figure(num="Polarization angle difference ratio")
-ax_dif_pa = fig_dif_pa.add_subplot(111, projection=wcs)
-cbar_ax_dif_pa = fig_dif_pa.add_axes([0.88, 0.12, 0.01, 0.75])
-ax_dif_pa.set_title(r"Degree of difference $\eta$ of the polarization from the 2 pipelines in the cut")
-im_dif_pa = ax_dif_pa.imshow(eta, vmin=0., vmax=2., origin='lower', cmap='bwr_r', label=r"$\eta$ between this pipeline and Kishimoto's")
-cbar_dif_pa = plt.colorbar(im_dif_pa, cax=cbar_ax_dif_pa,
-                           label=r"$\eta = \frac{2 \left|\theta_P^K-\theta_P^S\right|}{\sqrt{{\sigma^K_{\theta_P}}^2+{\sigma^S_{\theta_P}}^2}}$")
-ax_dif_pa.coords[0].set_axislabel('Right Ascension (J2000)')
-ax_dif_pa.coords[1].set_axislabel('Declination (J2000)')
-fig_dif_pa.savefig(path_join(root_dir_plot_S, "NGC1068_K_polang_diff.png"), bbox_inches="tight", dpi=300)
-
-# display both polarization maps to check consistency
-# plt.rcParams.update({'font.size': 15})
-fig = plt.figure(num="Polarization maps comparison", figsize=(10, 10))
-ax = fig.add_subplot(111, projection=wcs)
-fig.subplots_adjust(right=0.85)
-cbar_ax = fig.add_axes([0.88, 0.12, 0.01, 0.75])
-
-for d in [data_S, data_K]:
-    d['X'], d['Y'] = np.meshgrid(np.arange(d['I'].shape[1]), np.arange(d['I'].shape[0]))
-    d['xy_U'], d['xy_V'] = np.where(d['mask'], d['d_P']*np.cos(np.pi/2.+d['PA']), np.nan), np.where(d['mask'], d['d_P']*np.sin(np.pi/2.+d['PA']), np.nan)
-
-im0 = ax.imshow(data_S['I']*convert_flux, norm=LogNorm(data_S['I'][data_S['I'] > 0].min()*convert_flux, data_S['I']
-                [data_S['I'] > 0].max()*convert_flux), origin='lower', cmap='gray', label=r"$I_{STOKES}$ through this pipeline")
-quiv0 = ax.quiver(data_S['X'], data_S['Y'], data_S['xy_U'], data_S['xy_V'], units='xy', angles='uv', scale=0.5, scale_units='xy',
-                  pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.2, color='b', alpha=0.75, label="PA through this pipeline")
-quiv1 = ax.quiver(data_K['X'], data_K['Y'], data_K['xy_U'], data_K['xy_V'], units='xy', angles='uv', scale=0.5, scale_units='xy',
-                  pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.1, color='r', alpha=0.75, label="PA through Kishimoto's pipeline")
-
-ax.set_title(r"$SNR_P \geq$ "+str(SNRi_cut)+r"$\; & \; SNR_I \geq $"+str(SNRp_cut))
-# ax.coords.grid(True, color='white', ls='dotted', alpha=0.5)
-ax.coords[0].set_axislabel('Right Ascension (J2000)')
-ax.coords[0].set_axislabel_position('b')
-ax.coords[0].set_ticklabel_position('b')
-ax.coords[1].set_axislabel('Declination (J2000)')
-ax.coords[1].set_axislabel_position('l')
-ax.coords[1].set_ticklabel_position('l')
-# ax.axis('equal')
-
-cbar = plt.colorbar(im0, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
-ax.legend(loc='upper right')
-fig.savefig(path_join(root_dir_plot_S, "NGC1068_K_comparison.png"), bbox_inches="tight", dpi=300)
-
-# compute integrated polarization parameters on a specific cut
-for d in [data_S, data_K]:
-    d['I_dil'] = np.sum(d['I'][d['mask']])
-    d['sI_dil'] = np.sqrt(np.sum(d['sI'][d['mask']]**2))
-    d['Q_dil'] = np.sum(d['Q'][d['mask']])
-    d['sQ_dil'] = np.sqrt(np.sum(d['sQ'][d['mask']]**2))
-    d['U_dil'] = np.sum(d['U'][d['mask']])
-    d['sU_dil'] = np.sqrt(np.sum(d['sU'][d['mask']]**2))
-
-    d['P_dil'] = np.sqrt(d['Q_dil']**2+d['U_dil']**2)/d['I_dil']
-    d['sP_dil'] = np.sqrt((d['Q_dil']**2*d['sQ_dil']**2+d['U_dil']**2*d['sU_dil']**2)/(d['Q_dil']**2+d['U_dil']**2) +
-                          ((d['Q_dil']/d['I_dil'])**2+(d['U_dil']/d['I_dil'])**2)*d['sI_dil']**2)/d['I_dil']
-    d['d_P_dil'] = np.sqrt(d['P_dil']**2-d['sP_dil']**2)
-    d['PA_dil'] = princ_angle((90./np.pi)*np.arctan2(d['U_dil'], d['Q_dil']))
-    d['sPA_dil'] = princ_angle((90./(np.pi*(d['Q_dil']**2+d['U_dil']**2)))*np.sqrt(d['Q_dil']**2*d['sU_dil']**2+d['U_dil']**2*d['sU_dil']**2))
-print('From this pipeline :\n', "P = {0:.2f} ± {1:.2f} %\n".format(
-    data_S['d_P_dil']*100., data_S['sP_dil']*100.), "PA = {0:.2f} ± {1:.2f} °".format(data_S['PA_dil'], data_S['sPA_dil']))
-print("From Kishimoto's pipeline :\n", "P = {0:.2f} ± {1:.2f} %\n".format(
-    data_K['d_P_dil']*100., data_K['sP_dil']*100.), "PA = {0:.2f} ± {1:.2f} °".format(data_K['PA_dil'], data_K['sPA_dil']))
-
-# compare different types of error
-print("This pipeline : average sI/I={0:.2f} ; sQ/Q={1:.2f} ; sU/U={2:.2f} ; sP/P={3:.2f}".format(np.mean(data_S['sI'][data_S['mask']]/data_S['I'][data_S['mask']]), np.mean(
-    data_S['sQ'][data_S['mask']]/data_S['Q'][data_S['mask']]), np.mean(data_S['sU'][data_S['mask']]/data_S['U'][data_S['mask']]), np.mean(data_S['sP'][data_S['mask']]/data_S['P'][data_S['mask']])))
-print("Kishimoto's pipeline : average sI/I={0:.2f} ; sQ/Q={1:.2f} ; sU/U={2:.2f} ; sP/P={3:.2f}".format(np.mean(data_K['sI'][data_S['mask']]/data_K['I'][data_S['mask']]), np.mean(
-    data_K['sQ'][data_S['mask']]/data_K['Q'][data_S['mask']]), np.mean(data_K['sU'][data_S['mask']]/data_K['U'][data_S['mask']]), np.mean(data_K['sP'][data_S['mask']]/data_K['P'][data_S['mask']])))
-for d, i in zip(['I', 'Q', 'U', 'P', 'PA', 'sI', 'sQ', 'sU', 'sP', 'sPA'], [0, 1, 2, 5, 8, (3, 0, 0), (3, 1, 1), (3, 2, 2), 6, 9]):
-    data_K[d] = np.loadtxt(path_join(root_dir_K, d+'.txt'))
-    with fits.open(path_join(root_dir_data_S, filename_S)) as f:
-        if not type(i) is int:
-            data_S[d] = np.sqrt(f[i[0]].data[i[1], i[2]])
-        else:
-            data_S[d] = f[i].data
-    if i == 0:
-        header = f[i].header
-
-# from Kishimoto's pipeline : IQU_dir, IQU_shift, IQU_stat, IQU_trans
-# from my pipeline : raw_bg, raw_flat, raw_psf, raw_shift, raw_wav, IQU_dir
-#  but errors from my pipeline are propagated all along, how to compare then ?
-
-plt.show()

package/src/emission_center.py

@@ -0,0 +1,86 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+from pathlib import Path
+from sys import path as syspath
+
+syspath.append(str(Path(__file__).parent.parent))
+
+
+def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None):
+    from os.path import join as pathjoin
+
+    import numpy as np
+    from astropy.io.fits import open as fits_open
+    from astropy.wcs import WCS
+    from lib.plots import polarization_map
+    from lib.utils import CenterConf, PCconf
+    from matplotlib.patches import Rectangle
+    from matplotlib.pyplot import figure, show
+
+    output = []
+
+    Stokes = fits_open(infile)
+    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])],
+        [QN, UN, QN_ERR, UN_ERR],
+    ):
+        nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0]
+    Stokesconf = PCconf(QN, UN, QN_ERR, UN_ERR)
+    Stokesmask = Stokes["DATA_MASK"].data.astype(bool)
+    Stokessnr = np.zeros(Stokesmask.shape)
+    Stokessnr[Stokes["POL_DEG_ERR"].data > 0.0] = (
+        Stokes["POL_DEG_DEBIASED"].data[Stokes["POL_DEG_ERR"].data > 0.0] / Stokes["POL_DEG_ERR"].data[Stokes["POL_DEG_ERR"].data > 0.0]
+    )
+
+    if P_cut < 1.0:
+        Stokescentconf, Stokescenter = CenterConf(Stokesconf > P_cut, Stokes["POL_ANG"].data, Stokes["POL_ANG_ERR"].data)
+    else:
+        Stokescentconf, Stokescenter = CenterConf(Stokessnr > P_cut, Stokes["POL_ANG"].data, Stokes["POL_ANG_ERR"].data)
+    Stokespos = WCS(Stokes[0].header).pixel_to_world(*Stokescenter)
+
+    if target is 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)
+
+    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.contour(Stokescentconf, [0.01], colors="k", linewidths=3)
+    confcentcont = ax.contour(Stokescentconf, [0.01], colors="red")
+    # 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")
+    # 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")
+    # 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)
+
+    if output_dir is not None:
+        filename = pathjoin(output_dir, "%s_center.pdf" % target)
+        fig.savefig(filename, bbox_inches="tight", dpi=300, facecolor="None")
+        output.append(filename)
+    show()
+    return output
+
+
+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")
+    args = parser.parse_args()
+    exitcode = main(infile=args.file, P_cut=args.pcut, target=args.target, display=args.display, output_dir=args.output_dir)
+    print("Written to: ", exitcode)

package/src/get_cdelt.py

@@ -1,4 +1,9 @@
-#!/usr/bin/python
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+from pathlib import Path
+from sys import path as syspath
+
+syspath.append(str(Path(__file__).parent.parent))
 
 
 def main(infiles=None):
@@ -13,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

@@ -0,0 +1,56 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+from pathlib import Path
+from sys import path as syspath
+
+syspath.append(str(Path(__file__).parent.parent))
+
+import numpy as np
+from astropy.io import fits
+from lib.plots import overplot_pol, overplot_radio
+from matplotlib.colors import LogNorm
+
+Stokes_UV = fits.open("./data/IC5063/5918/IC5063_FOC_b0.10arcsec_c0.20arcsec.fits")
+# Stokes_18GHz = fits.open("./data/IC5063/radio/IC5063_18GHz.fits")
+# Stokes_24GHz = fits.open("./data/IC5063/radio/IC5063_24GHz.fits")
+# Stokes_103GHz = fits.open("./data/IC5063/radio/IC5063_103GHz.fits")
+# Stokes_229GHz = fits.open("./data/IC5063/radio/IC5063_229GHz.fits")
+# Stokes_357GHz = fits.open("./data/IC5063/radio/IC5063_357GHz.fits")
+# Stokes_S2 = fits.open("./data/IC5063/POLARIZATION_COMPARISON/S2_rot_crop.fits")
+Stokes_IR = fits.open("./data/IC5063/IR/u2e65g01t_c0f_rot.fits")
+
+# levelsMorganti = np.array([1.,2.,3.,8.,16.,32.,64.,128.])
+# levelsMorganti = np.logspace(-0.1249, 1.97, 7) / 100.0
+
+# levels18GHz = levelsMorganti * Stokes_18GHz[0].data.max()
+# A = overplot_radio(Stokes_UV, Stokes_18GHz)
+# A.plot(levels=levels18GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/18GHz_overplot.pdf", scale_vec=None)
+
+# levels24GHz = levelsMorganti * Stokes_24GHz[0].data.max()
+# B = overplot_radio(Stokes_UV, Stokes_24GHz)
+# B.plot(levels=levels24GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/24GHz_overplot.pdf", scale_vec=None)
+
+# levels103GHz = levelsMorganti * Stokes_103GHz[0].data.max()
+# C = overplot_radio(Stokes_UV, Stokes_103GHz)
+# C.plot(levels=levels103GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/103GHz_overplot.pdf", scale_vec=None)
+
+# levels229GHz = levelsMorganti * Stokes_229GHz[0].data.max()
+# D = overplot_radio(Stokes_UV, Stokes_229GHz)
+# D.plot(levels=levels229GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/229GHz_overplot.pdf", scale_vec=None)
+
+# levels357GHz = levelsMorganti * Stokes_357GHz[0].data.max()
+# E = overplot_radio(Stokes_UV, Stokes_357GHz)
+# E.plot(levels=levels357GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/357GHz_overplot.pdf", scale_vec=None)
+
+# F = overplot_pol(Stokes_UV, Stokes_S2)
+# F.plot(P_cut=3.0, SNRi_cut=80.0, savename='./plots/IC5063/S2_overplot.pdf', norm=LogNorm(vmin=5e-20,vmax=5e-18))
+
+G = overplot_pol(Stokes_UV, Stokes_IR, cmap="inferno")
+G.plot(
+    P_cut=0.99,
+    SNRi_cut=1.0,
+    savename="./plots/IC5063/IR_overplot.pdf",
+    scale_vec=None,
+    norm=LogNorm(Stokes_IR[0].data.max() * Stokes_IR[0].header["photflam"] / 1e3, Stokes_IR[0].data.max() * Stokes_IR[0].header["photflam"]),
+    cmap="inferno",
+)

package/src/overplot_MRK463E.py

@@ -0,0 +1,46 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+from pathlib import Path
+from sys import path as syspath
+
+syspath.append(str(Path(__file__).parent.parent))
+
+import numpy as np
+from astropy.io import fits
+from lib.plots import overplot_chandra, overplot_pol
+from matplotlib.colors import LogNorm
+
+Stokes_UV = fits.open("./data/MRK463E/5960/MRK463E_FOC_b0.05arcsec_c0.10arcsec.fits")
+# Stokes_IR = fits.open("./data/MRK463E/WFPC2/IR_rot_crop.fits")
+# Stokes_Xr = fits.open("./data/MRK463E/Chandra/X_ray_crop.fits")
+Radio = fits.open("./data/MRK463E/EMERLIN/Voorwerpjes_1356+1822_1356+1822_uniform-image.fits")
+
+# levels = np.geomspace(1.0, 99.0, 7)
+
+# A = overplot_chandra(Stokes_UV, Stokes_Xr, norm=LogNorm())
+# A.plot(levels=levels, P_cut=0.99, SNRi_cut=1.0, scale_vec=5, zoom=1, savename="./plots/MRK463E/Chandra_overplot.pdf")
+# A.write_to(path1="./data/MRK463E/FOC_data_Chandra.fits", path2="./data/MRK463E/Chandra_data.fits", suffix="aligned")
+
+# levels = np.array([0.8, 2, 5, 10, 20, 50]) / 100.0 * Stokes_UV[0].header["photflam"]
+# B = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm())
+# 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([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.99,
+    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,
+    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

@@ -0,0 +1,32 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+from pathlib import Path
+from sys import path as syspath
+
+syspath.append(str(Path(__file__).parent.parent))
+
+import numpy as np
+from astropy.io import fits
+from lib.plots import overplot_pol, plt
+from matplotlib.colors import LogNorm
+
+Stokes_UV = fits.open("./data/NGC1068/5144/NGC1068_FOC_b0.05arcsec_c0.07arcsec_crop.fits")
+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.add_vector(
+    A.other_wcs.celestial.wcs.crpix - (1.0, 1.0),
+    pol_deg=0.124,
+    pol_ang=100.7,
+    width=2.0,
+    linewidth=1.0,
+    scale=1.0 / (A.px_scale * 6.0),
+    edgecolor="w",
+    color="b",
+    label=r"IXPE torus: P = $12.4 (\pm 3.6)$%, PA = $100.7 (\pm 8.3)$°",
+)
+A.fig_overplot.savefig("./plots/NGC1068/NGC1068_radio_overplot_full.pdf", dpi=300, bbox_inches="tight")
+plt.show()
+A.write_to(path1="./data/NGC1068/FOC_Radio.fits", path2="./data/NGC1068/Radio_FOC.fits", suffix="aligned")

package/src/readfos.py

@@ -0,0 +1,850 @@
+#!/usr/bin/env python3
+import matplotlib.pyplot as plt
+import numpy as np
+from astropy.io.fits import getheader, getdata, hdu
+from os.path import join as join_path, exists as path_exists
+from os import system
+from copy import deepcopy
+
+# consecutive spectra are made up of the summ of all previous ACCUMs, so the S/N increases along sequence
+#  _c0f.fits - calibrated vacuum wavelength
+#  _c1f.fits - calibrated fluxes (ergs sec^-1 cm^-2 Angs^-1)
+#  _c2f.fits - statistical errors (no sky, bkg subtraction, flatfield or sensitivity error)
+#  _c3f.fits - for SPECTROPOLARIMETRY mode contains I, Q, U, V, linear and circular polarization and polarization position angle
+#  _c4f.fits - object+sky count rate spectrum (corrected for overscanning, noise rejection, lost signal)
+#  _c5f.fits - flat-fielded object count rate spectrum (corrected for paired pulses, detector background, flatfield structure, GIM effects)
+#  _c6f.fits - flat-fielded sky count rate spectrum (corrected for paired pulses, detector background, flatfield structure, GIM effects)
+#  _c7f.fits - background count rate spectrum (scaled background subtracted from c4 products)
+#  _c8f.fits - flat-fielded and sky-subtracted object count rate spectrum
+
+
+def princ_angle(ang):
+    """
+    Return the principal angle in the 0° to 180° quadrant as PA is always
+    defined at p/m 180°.
+    ----------
+    Inputs:
+    ang : float, numpy.ndarray
+        Angle in degrees. Can be an array of angles.
+    ----------
+    Returns:
+    princ_ang : float, numpy.ndarray
+        Principal angle in the 0°-180° quadrant in the same shape as input.
+    """
+    if not isinstance(ang, np.ndarray):
+        A = np.array([ang])
+    else:
+        A = np.array(ang)
+    while np.any(A < 0.0):
+        A[A < 0.0] = A[A < 0.0] + 360.0
+    while np.any(A >= 180.0):
+        A[A >= 180.0] = A[A >= 180.0] - 180.0
+    if type(ang) is type(A):
+        return A
+    else:
+        return A[0]
+
+
+class specpol(object):
+    """
+    Class object for studying spectropolarimetry.
+    """
+
+    def __init__(self, other=None):
+        if isinstance(other, __class__):
+            # Copy constructor
+            self.hd = deepcopy(other.hd)
+            self.bin_edges = deepcopy(other.bin_edges)
+            self.wav = deepcopy(other.wav)
+            self.wav_err = deepcopy(other.wav_err)
+            self.I = deepcopy(other.I)
+            self.Q = deepcopy(other.Q)
+            self.U = deepcopy(other.U)
+            self.V = deepcopy(other.V)
+            self.IQUV_cov = deepcopy(other.IQUV_cov)
+            if hasattr(other, "I_r"):
+                self.I_r = deepcopy(other.I_r)
+                self.I_r_err = deepcopy(other.I_r_err)
+                self.wav_r = deepcopy(other.wav_r)
+                self.wav_r_err = deepcopy(other.wav_r_err)
+        elif isinstance(other, str):
+            self.from_txt(other)
+        else:
+            # Initialise to zero
+            if isinstance(other, int):
+                self.zero(other)
+            else:
+                self.zero()
+
+    @classmethod
+    def zero(self, n=1):
+        """
+        Set all values to zero.
+        """
+        self.hd = dict([])
+        self.hd["TARGNAME"], self.hd["PROPOSID"], self.hd["ROOTNAME"], self.hd["APER_ID"] = "", 0, "", ""
+        self.hd["DENSITY"] = False
+        self.hd["XUNIT"], self.hd["YUNIT"] = r"Wavelength [m]", r"{0:s}F [$10^{{{1:d}}}$ count s$^{{-1}}$]"
+        self.bin_edges = np.zeros(n + 1)
+        self.wav = np.zeros(n)
+        self.wav_err = np.zeros((n, 2))
+        self.I = np.zeros(n)
+        self.Q = np.zeros(n)
+        self.U = np.zeros(n)
+        self.V = np.zeros(n)
+        self.IQUV_cov = np.zeros((4, 4, n))
+
+    def rest(self, wav=None, z=None):
+        if z is None and self.hd["TARGNAME"] == "":
+            z = 0
+        elif z is None and "REDSHIFT" not in self.hd.keys():
+            from astroquery.ipac.ned import Ned
+
+            z = Ned.query_object(self.hd["TARGNAME"])["Redshift"][0]
+            self.hd["REDSHIFT"] = z
+        elif z is None:
+            z = self.hd["REDSHIFT"]
+        if wav is None:
+            wav = self.wav
+        return wav / (z + 1)
+
+    def unrest(self, wav=None, z=None):
+        if z is None and self.hd["TARGNAME"] == "":
+            z = 0
+        elif z is None and "REDSHIFT" not in self.hd.keys():
+            from astroquery.ipac.ned import Ned
+
+            z = Ned.query_object(self.hd["TARGNAME"])["Redshift"][0]
+            self.hd["REDSHIFT"] = z
+        elif z is None:
+            z = self.hd["REDSHIFT"]
+        if wav is None:
+            wav = self.wav
+        return wav * (z + 1)
+
+    @property
+    def I_err(self):
+        return np.sqrt(self.IQUV_cov[0][0])
+
+    @property
+    def Q_err(self):
+        return np.sqrt(self.IQUV_cov[1][1])
+
+    @property
+    def U_err(self):
+        return np.sqrt(self.IQUV_cov[2][2])
+
+    @property
+    def V_err(self):
+        return np.sqrt(self.IQUV_cov[3][3])
+
+    @property
+    def QN(self):
+        np.seterr(divide="ignore", invalid="ignore")
+        return self.Q / np.where(self.I > 0, self.I, np.nan)
+
+    @property
+    def QN_err(self):
+        np.seterr(divide="ignore", invalid="ignore")
+        return self.Q_err / np.where(self.I > 0, self.I, np.nan)
+
+    @property
+    def UN(self):
+        np.seterr(divide="ignore", invalid="ignore")
+        return self.U / np.where(self.I > 0, self.I, np.nan)
+
+    @property
+    def UN_err(self):
+        np.seterr(divide="ignore", invalid="ignore")
+        return self.U_err / np.where(self.I > 0, self.I, np.nan)
+
+    @property
+    def VN(self):
+        np.seterr(divide="ignore", invalid="ignore")
+        return self.V / np.where(self.I > 0, self.I, np.nan)
+
+    @property
+    def VN_err(self):
+        np.seterr(divide="ignore", invalid="ignore")
+        return self.V_err / np.where(self.I > 0, self.I, np.nan)
+
+    @property
+    def PF(self):
+        np.seterr(divide="ignore", invalid="ignore")
+        return np.sqrt(self.Q**2 + self.U**2 + self.V**2)
+
+    @property
+    def PF_err(self):
+        np.seterr(divide="ignore", invalid="ignore")
+        return np.sqrt(self.Q**2 * self.Q_err**2 + self.U**2 * self.U_err**2 + self.V**2 * self.V_err**2) / np.where(self.PF > 0, self.PF, np.nan)
+
+    @property
+    def P(self):
+        np.seterr(divide="ignore", invalid="ignore")
+        return self.PF / np.where(self.I > 0, self.I, np.nan)
+
+    @property
+    def P_err(self):
+        np.seterr(divide="ignore", invalid="ignore")
+        return np.where(self.I > 0, np.sqrt(self.PF_err**2 + (self.PF / self.I) ** 2 * self.I_err**2) / self.I, np.nan)
+
+    @property
+    def PA(self):
+        return princ_angle((90.0 / np.pi) * np.arctan2(self.U, self.Q))
+
+    @property
+    def PA_err(self):
+        return princ_angle((90.0 / np.pi) * np.sqrt(self.U**2 * self.Q_err**2 + self.Q**2 * self.U_err**2) / np.where(self.PF > 0, self.PF**2, np.nan))
+
+    def rotate(self, angle):
+        alpha = np.pi / 180.0 * angle
+        mrot = np.array(
+            [
+                [1.0, 0.0, 0.0, 0.0],
+                [0.0, np.cos(2.0 * alpha), np.sin(2.0 * alpha), 0.0],
+                [0.0, -np.sin(2.0 * alpha), np.cos(2.0 * alpha), 0.0],
+                [0.0, 0.0, 0.0, 1.0],
+            ]
+        )
+        self.I, self.Q, self.U, self.V = np.dot(mrot, np.array([self.I, self.Q, self.U, self.V]))
+        self.IQUV_cov = np.dot(mrot, np.dot(self.IQUV_cov.T, mrot.T).T)
+
+    def bin(self, bin_edges):
+        """
+        Rebin spectra to given list of bin edges.
+        """
+        # Get new binning distribution and define new empty spectra
+        in_bin = np.digitize(self.wav, bin_edges) - 1
+        out = specpol(bin_edges.shape[0] - 1)
+        if hasattr(self, "I_r"):
+            # Propagate "raw" flux spectra to new bin
+            out.I_r = deepcopy(self.I_r)
+            out.I_r_err = deepcopy(self.I_r_err)
+            out.wav_r = deepcopy(self.wav_r)
+            out.wav_r_err = deepcopy(self.wav_r_err)
+        else:
+            # Create "raw" flux spectra from previously unbinned spectra
+            out.I_r = deepcopy(self.I[self.I > 0.0])
+            out.I_r_err = deepcopy(self.I_err[self.I > 0.0])
+            out.wav_r = deepcopy(self.wav[self.I > 0.0])
+            out.wav_r_err = deepcopy(self.wav_err[self.I > 0.0])
+
+        for i in range(bin_edges.shape[0] - 1):
+            # Set the wavelength as the mean wavelength of acquisitions in bin, default to the bin center
+            out.wav[i] = np.mean(self.wav[in_bin == i]) if np.any(in_bin == i) else 0.5 * (bin_edges[i] + bin_edges[i + 1])
+            out.wav_err[i] = (out.wav[i] - bin_edges[i], bin_edges[i + 1] - out.wav[i])
+
+            if self.hd["DENSITY"] and np.any(in_bin == i):
+                # If flux density, convert to flux before converting back to the new density
+                wav1 = np.abs(self.wav_err[in_bin == i]).sum(axis=1)
+                wav2 = np.abs(out.wav_err[i]).sum()
+            else:
+                wav1, wav2 = 1.0, 1.0
+            out.I[i] = np.sum(self.I[in_bin == i] * wav1) / wav2 if np.any(in_bin == i) else 0.0
+            out.Q[i] = np.sum(self.Q[in_bin == i] * wav1) / wav2 if np.any(in_bin == i) else 0.0
+            out.U[i] = np.sum(self.U[in_bin == i] * wav1) / wav2 if np.any(in_bin == i) else 0.0
+            out.V[i] = np.sum(self.V[in_bin == i] * wav1) / wav2 if np.any(in_bin == i) else 0.0
+            for m in range(4):
+                # Quadratically sum the uncertainties
+                out.IQUV_cov[m][m][i] = np.sum(self.IQUV_cov[m][m][in_bin == i] * wav1**2) / wav2**2 if np.any(in_bin == i) else 0.0
+                for n in [k for k in range(4) if k != m]:
+                    out.IQUV_cov[m][n][i] = np.sqrt(np.sum((self.IQUV_cov[m][n][in_bin == i] * wav1) ** 2)) / wav2 if np.any(in_bin == i) else 0.0
+        # Update bin edges and header
+        out.bin_edges = bin_edges
+        out.hd = deepcopy(self.hd)
+        out.hd["NAXIS1"] = bin_edges.shape[0] - 1
+        out.hd["DATAMIN"], out.hd["DATAMAX"] = out.I.min(), out.I.max()
+        out.hd["MINWAV"], out.hd["MAXWAV"] = out.wav.min(), out.wav.max()
+        out.hd["STEPWAV"] = np.max(bin_edges[1:] - bin_edges[:-1])
+        return out
+
+    def bin_size(self, size):
+        """
+        Rebin spectra to selected bin size in Angstrom.
+        """
+        bin_edges = np.arange(self.bin_edges.min(), self.bin_edges.max() + size, size, dtype=np.float32)
+        return self.bin(bin_edges)
+
+    def from_txt(self, filename, data_dir=""):
+        """
+        Fill current spectra from a text file.
+        """
+        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_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):
+            self.IQUV_cov[i][i] = deepcopy(self.IQUV_cov[i][i]) ** 2
+        with open(join_path(data_dir, filename)) as f:
+            self.hd["TARGNAME"], self.hd["PROPOSID"], self.hd["ROOTNAME"], self.hd["APER_ID"], self.hd["XUNIT"], self.hd["YUNIT"] = f.readline()[2:].split(";")
+
+    def dump_txt(self, filename, output_dir=""):
+        """
+        Dump current spectra to a text file.
+        """
+        header = ";".join([self.hd["TARGNAME"], str(self.hd["PROPOSID"]), self.hd["ROOTNAME"], self.hd["APER_ID"], self.hd["XUNIT"], self.hd["YUNIT"]])
+        header += "\nwav\t wav_err\t I\t I_err\t Q\t Q_err\t U\t U_err\t V\t V_err\t P\t P_err\t PA\t PA_err"
+        data_dump = np.array(
+            [
+                self.wav,
+                self.wav_err.mean(axis=1),
+                self.I,
+                self.I_err,
+                self.Q,
+                self.Q_err,
+                self.U,
+                self.U_err,
+                self.V,
+                self.V_err,
+                self.P,
+                self.P_err,
+                self.PA,
+                self.PA_err,
+            ]
+        ).T
+        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=""):
+        """
+        Display current spectra.
+        """
+        if fig is None:
+            plt.rcParams.update({"font.size": 15})
+            if ax is None:
+                self.fig, self.ax = plt.subplots(3, 1, sharex=True, figsize=(20, 15))
+                self.fig.subplots_adjust(hspace=0)
+                self.fig.suptitle("_".join([self.hd["TARGNAME"], str(self.hd["PROPOSID"]), self.hd["ROOTNAME"], self.hd["APER_ID"]]))
+            else:
+                self.ax = ax
+        else:
+            if ax is None:
+                self.fig = fig
+                self.ax = self.fig.add_subplot(111)
+            else:
+                self.fig = fig
+                self.ax = ax
+        if isinstance(self.ax, np.ndarray):
+            if self.ax.shape[0] == 2:
+                ax1, ax2 = self.ax[:2]
+                ax22 = ax2.twinx()
+                ax2.set_xlabel(self.hd["XUNIT"])
+                secax1 = ax1.secondary_xaxis("top", functions=(self.rest, self.unrest))
+                secax1.set_xlabel(r"Rest " + self.hd["XUNIT"])
+            else:
+                ax1, ax2, ax22 = self.ax[::-1]
+        else:
+            ax1 = self.ax
+
+        # Display flux and polarized flux on first ax
+        if hasattr(self, "I_r"):
+            # If available, display "raw" total flux
+            yoffset = np.floor(np.log10(self.I_r[self.I_r > 0.0].min())).astype(int)
+            yoff = 10.0**yoffset
+            ymin, ymax = (
+                np.min((self.I_r - 1.5 * self.I_r_err)[self.I_r > 1.5 * self.I_r_err]) / yoff,
+                np.max((self.I_r + self.I_r_err * 1.5)[self.I_r > 1.5 * self.I_r_err]) / yoff,
+            )
+            xmin, xmax = np.min(self.wav_r - self.wav_r_err[:, 0]), np.max(self.wav_r + self.wav_r_err[:, 1])
+            ax1.errorbar(self.wav_r, self.I_r / yoff, xerr=self.wav_r_err.T, yerr=self.I_r_err / yoff, color="k", fmt=".", label="I")
+        else:
+            yoffset = np.floor(np.log10(self.I[self.I > 0.0].min())).astype(int)
+            yoff = 10.0**yoffset
+            ymin, ymax = (
+                np.min((self.I - 1.5 * self.I_err)[self.I > 1.5 * self.I_err]) / yoff,
+                np.max((self.I + self.I_err * 1.5)[self.I > 1.5 * self.I_err]) / yoff,
+            )
+            xmin, xmax = np.min(self.wav - self.wav_err[:, 0]), np.max(self.wav + self.wav_err[:, 1])
+            ax1.errorbar(self.wav, self.I / yoff, xerr=self.wav_err.T, yerr=self.I_err / yoff, color="k", fmt=".", label="I")
+
+        ax1.set_xlim([np.min([xmin, self.bin_edges.min()]), np.max([xmax, self.bin_edges.max()])])
+        ax1.set_xlabel(self.hd["XUNIT"])
+        ax1.set_ylim([ymin, ymax])
+        ax1.set_ylabel(self.hd["YUNIT"].format("", yoffset))
+
+        # ax11 = ax1.twinx()
+        # pfoffset = np.floor(np.log10(self.PF[self.PF > 0.0].min())).astype(int)
+        # pfoff = 10.0**pfoffset
+        # ax11.errorbar(self.wav, self.PF / pfoff, xerr=self.wav_err.T, yerr=self.PF_err / pfoff, color="b", fmt=".", label="PF")
+        # ax11.set_ylim(
+        #     [
+        #         ymin * yoff * self.P[np.logical_and(self.P > 0.0, np.isfinite(self.P))].min() / pfoff,
+        #         ymax * yoff * self.P[np.logical_and(self.P > 0.0, np.isfinite(self.P))].max() / pfoff,
+        #     ]
+        # )
+        # ax11.set_ylabel(self.hd["YUNIT"].format("Px", pfoffset), color="b")
+        # ax11.tick_params(axis="y", color="b", labelcolor="b")
+
+        # ax1.legend(ncols=2, loc=1)
+
+        if isinstance(self.ax, np.ndarray):
+            # When given 2 axes, display P and PA on second
+            ax2.errorbar(self.wav, self.P * 100.0, xerr=self.wav_err.T, yerr=self.P_err * 100.0, color="b", fmt=".", label="P")
+            pmin, pmax = (
+                np.min(self.P[self.I > 0.0] - 1.5 * self.P_err[self.I > 0.0]) * 100.0,
+                np.max(self.P[self.I > 0.0] + 1.5 * self.P_err[self.I > 0.0]) * 100.0,
+            )
+            ax2.set_ylim([pmin if pmin > 0.0 else 0.0, pmax if pmax < 100.0 else 100.0])
+            ax2.set_ylabel(r"P [%]", color="b")
+            ax2.tick_params(axis="y", color="b", labelcolor="b")
+
+            ax22.errorbar(self.wav, self.PA, xerr=self.wav_err.T, yerr=self.PA_err, color="r", fmt=".", label="PA [°]")
+            pamin, pamax = np.min(self.PA[self.I > 0.0] - 1.5 * self.PA_err[self.I > 0.0]), np.max(self.PA[self.I > 0.0] + 1.5 * self.PA_err[self.I > 0.0])
+            ax22.set_ylim([pamin if pamin > 0.0 else 0.0, pamax if pamax < 180.0 else 180.0])
+            ax22.set_ylabel(r"PA [°]", color="r")
+            ax22.tick_params(axis="y", color="r", labelcolor="r")
+
+            secax22 = ax22.secondary_xaxis("top", functions=(self.rest, self.unrest))
+            secax22.set_xlabel(r"Rest " + self.hd["XUNIT"])
+            h2, l2 = ax2.get_legend_handles_labels()
+            h22, l22 = ax22.get_legend_handles_labels()
+            if self.ax.shape[0] == 2:
+                ax2.legend(h2 + h22, l2 + l22, ncols=2, loc=1)
+
+        if hasattr(self, "fig") and savename is not None:
+            self.fig.savefig(join_path(plots_folder, savename + ".pdf"), dpi=300, bbox_inches="tight")
+            return self.fig, self.ax, join_path(plots_folder, savename + ".pdf")
+        elif hasattr(self, "fig"):
+            return self.fig, self.ax
+        else:
+            return self.ax
+
+    def __add__(self, other):
+        """
+        Spectra addition, if not same binning concatenate both spectra binning.
+        """
+        if (self.bin_edges.shape == other.bin_edges.shape) and np.all(self.bin_edges == other.bin_edges):
+            bin_edges = deepcopy(self.bin_edges)
+        else:
+            # If different binning, concatenate binnings
+            if self.bin_edges[0] <= other.bin_edges[0]:
+                bin_edges = deepcopy(self.bin_edges)
+            else:
+                bin_edges = deepcopy(other.bin_edges)
+            if other.bin_edges[-1] > bin_edges[-1]:
+                bin_edges = np.concat((bin_edges, deepcopy(other.bin_edges[other.bin_edges > bin_edges[-1]])), axis=0)
+            elif self.bin_edges[-1] > bin_edges[-1]:
+                bin_edges = np.concat((bin_edges, deepcopy(self.bin_edges[self.bin_edges > bin_edges[-1]])), axis=0)
+        # Rebin spectra to be added to ensure same binning
+        spec_a = specpol(specpol(self).bin(bin_edges=bin_edges))
+        spec_b = specpol(specpol(other).bin(bin_edges=bin_edges))
+
+        # Create sum spectra
+        spec = specpol(bin_edges.shape[0] - 1)
+        spec.hd = deepcopy(self.hd)
+        spec.bin_edges = bin_edges
+        spec.wav = np.mean([spec_a.wav, spec_b.wav], axis=0)
+        spec.wav_err = np.array([spec.wav - spec.bin_edges[:-1], spec.bin_edges[1:] - spec.wav]).T
+
+        # Propagate "raw" flux spectra to sum
+        if hasattr(self, "I_r") and hasattr(other, "I_r"):
+            # Deal with the concatenation of the "raw" total flux spectra
+            if self.wav_r[0] <= other.wav_r[0]:
+                inter = other.wav_r[0], self.wav_r[-1]
+                spec.wav_r = deepcopy(np.concat((self.wav_r, other.wav_r[other.wav_r > self.wav_r[-1]])))
+                spec.wav_r_err = deepcopy(np.concat((self.wav_r_err, other.wav_r_err[other.wav_r > self.wav_r[-1]]), axis=0))
+                spec.I_r = deepcopy(np.concat((self.I_r, other.I_r[other.wav_r > self.wav_r[-1]])))
+                spec.I_r_err = deepcopy(np.concat((self.I_r_err, other.I_r_err[other.wav_r > self.wav_r[-1]]), axis=0))
+            else:
+                inter = self.wav_r[0], other.wav_r[-1]
+                spec.wav_r = deepcopy(np.concat((other.wav_r, self.wav_r[self.wav_r > other.wav_r[-1]])))
+                spec.wav_r_err = deepcopy(np.concat((other.wav_r_err, self.wav_r_err[self.wav_r > other.wav_r[-1]]), axis=0))
+                spec.I_r = deepcopy(np.concat((other.I_r, self.I_r[self.wav_r > other.wav_r[-1]])))
+                spec.I_r_err = deepcopy(np.concat((other.I_r_err, self.I_r_err[self.wav_r > other.wav_r[-1]]), axis=0))
+            # When both spectra intersect, compute intersection as the mean
+            edges = np.concat((spec.wav_r - spec.wav_r_err[:, 0], [spec.wav_r[-1] + spec.wav_r_err[-1, 1]]))
+            edges.sort()
+            bin, bino = np.digitize(self.wav_r, edges) - 1, np.digitize(other.wav_r, edges) - 1
+            for w in np.arange(spec.wav_r.shape[0])[np.logical_and(spec.wav_r >= inter[0], spec.wav_r <= inter[1])]:
+                if self.hd["DENSITY"] and np.any(bin == w):
+                    # If flux density, convert to flux before converting back to the new density
+                    wav, wavo = (
+                        np.abs(self.wav_r_err[bin == w]).sum(axis=1) * (self.I_r[bin == w] > self.I_r_err[bin == w]),
+                        np.abs(other.wav_r_err[bino == w]).sum(axis=1) * (other.I_r[bino == w] > other.I_r_err[bino == w]),
+                    )
+                    wavs = np.abs(spec.wav_r_err[w]).sum()
+                else:
+                    wav, wavo, wavs = 1.0, 1.0, 1.0
+                n = np.sum(self.I_r[bin == w] > self.I_r_err[bin == w]) + np.sum(other.I_r[bino == w] > other.I_r_err[bino == w])
+                spec.I_r[w] = np.sum(np.concat([self.I_r[bin == w] * wav, other.I_r[bino == w] * wavo])) / wavs / n
+                spec.I_r_err[w] = np.sqrt(np.sum(np.concat([self.I_r_err[bin == w] ** 2 * wav**2, other.I_r_err[bino == w] ** 2 * wavo**2]))) / wavs / n
+
+        # Sum stokes fluxes
+        spec.I = deepcopy(spec_a.I + spec_b.I)
+        spec.Q = deepcopy(spec_a.Q + spec_b.Q)
+        spec.U = deepcopy(spec_a.U + spec_b.U)
+        spec.V = deepcopy(spec_a.V + spec_b.V)
+        # Quadratically sum uncertainties
+        for i in range(4):
+            spec.IQUV_cov[i][i] = deepcopy(spec_a.IQUV_cov[i][i] + spec_b.IQUV_cov[i][i])
+            for j in [k for k in range(4) if k != i]:
+                spec.IQUV_cov[i][j] = deepcopy(np.sqrt(spec_a.IQUV_cov[i][j] ** 2 + spec_b.IQUV_cov[i][j] ** 2))
+
+        # Update header to reflect sum
+        spec.hd["DATAMIN"], spec.hd["DATAMAX"] = spec.I.min(), spec.I.max()
+        spec.hd["MINWAV"], spec.hd["MAXWAV"] = spec.wav.min(), spec.wav.max()
+        spec.hd["EXPTIME"] = spec_a.hd["EXPTIME"] + spec_b.hd["EXPTIME"]
+        rootnames = [spec_a.hd["ROOTNAME"], spec_b.hd["ROOTNAME"]]
+        spec.hd["ROOTNAME"] = "".join(p for p, *r in zip(*rootnames) if all(p == c for c in r)) + "_SUM"
+        return spec
+
+    def __deepcopy__(self, memo={}):
+        spec = specpol(self)
+        spec.__dict__.update(self.__dict__)
+
+        spec.hd = deepcopy(self.hd, memo)
+        spec.bin_edges = deepcopy(spec.bin_edges, memo)
+        spec.wav = deepcopy(self.wav, memo)
+        spec.wav_err = deepcopy(self.wav_err, memo)
+        spec.I = deepcopy(self.I, memo)
+        spec.Q = deepcopy(self.Q, memo)
+        spec.U = deepcopy(self.U, memo)
+        spec.V = deepcopy(self.V, memo)
+        spec.IQUV_cov = deepcopy(self.IQUV_cov, memo)
+
+        return spec
+
+
+class FOSspecpol(specpol):
+    """
+    Class object for studying FOS SPECTROPOLARYMETRY mode spectra.
+    """
+
+    def __init__(self, stokes, data_folder=""):
+        """
+        Initialise object from fits filename, fits hdulist or copy.
+        """
+        if isinstance(stokes, __class__):
+            # Copy constructor
+            self.rootname = deepcopy(stokes.rootname)
+            self.hd = deepcopy(stokes.hd)
+            self.bin_edges = deepcopy(stokes.bin_edges)
+            self.wav = deepcopy(stokes.wav)
+            self.wav_err = deepcopy(stokes.wav_err)
+            self.I = deepcopy(stokes.I)
+            self.Q = deepcopy(stokes.Q)
+            self.U = deepcopy(stokes.U)
+            self.V = deepcopy(stokes.V)
+            self.IQUV_cov = deepcopy(stokes.IQUV_cov)
+            self.P_fos = deepcopy(stokes.P_fos)
+            self.P_fos_err = deepcopy(stokes.P_fos_err)
+            self.PC_fos = deepcopy(stokes.PC_fos)
+            self.PC_fos_err = deepcopy(stokes.PC_fos_err)
+            self.PA_fos = deepcopy(stokes.PA_fos)
+            self.PA_fos_err = deepcopy(stokes.PA_fos_err)
+            self.subspec = {}
+            for name in ["PASS1", "PASS2", "PASS12", "PASS12corr"]:
+                spec = deepcopy(stokes.subspec[name])
+                self.subspec[name] = spec
+        elif stokes is None or isinstance(stokes, int):
+            self.zero(n=stokes)
+        else:
+            self.from_file(stokes, data_folder=data_folder)
+
+    @classmethod
+    def zero(self, n=1):
+        """
+        Set all values to zero.
+        """
+        self.rootname = ""
+        self.hd = dict([])
+        self.hd["DENSITY"] = True
+        self.hd["TARGNAME"] = "Undefined"
+        self.hd["XUNIT"], self.hd["YUNIT"] = r"Wavelength [$\AA$]", r"{0:s}F$_\lambda$ [$10^{{{1:d}}}$ erg s$^{{-1}}$ cm$^{{-2}} \AA^{{-1}}$]"
+        self.bin_edges = np.zeros((4, n + 1))
+        self.wav = np.zeros((4, n))
+        self.wav_err = np.zeros((4, n, 2))
+        self.I = np.zeros((4, n))
+        self.Q = np.zeros((4, n))
+        self.U = np.zeros((4, n))
+        self.V = np.zeros((4, n))
+        self.IQUV_cov = np.zeros((4, 4, 4, n))
+
+        self.subspec = {}
+        for i, name in enumerate(["PASS1", "PASS2", "PASS12", "PASS12corr"]):
+            spec = specpol(n)
+            spec.hd, spec.wav, spec.wav_err, spec.I, spec.Q, spec.U, spec.V = self.hd, self.wav[i], self.wav_err[i], self.I[i], self.Q[i], self.U[i], self.V[i]
+            spec.IQUV_cov = self.IQUV_cov[:, :, i, :]
+            self.subspec[name] = spec
+
+        self.P_fos = np.zeros(self.I.shape)
+        self.P_fos_err = np.zeros(self.I.shape)
+        self.PC_fos = np.zeros(self.I.shape)
+        self.PC_fos_err = np.zeros(self.I.shape)
+        self.PA_fos = np.zeros(self.I.shape)
+        self.PA_fos_err = np.zeros(self.I.shape)
+
+    def from_file(self, stokes, data_folder=""):
+        """
+        Initialise object from fits file or hdulist.
+        """
+        if isinstance(stokes, str):
+            self.rootname = stokes.split("_")[0]
+            self.hd = dict(getheader(join_path(data_folder, self.rootname + "_c3f.fits")))
+            wav = getdata(join_path(data_folder, self.rootname + "_c0f.fits"))
+            stokes = getdata(join_path(data_folder, self.rootname + "_c3f.fits"))
+        elif isinstance(stokes, hdu.hdulist.HDUList):
+            self.hd = dict(stokes.header)
+            self.rootname = self.hd["FILENAME"].split("_")[0]
+            wav = getdata(join_path(data_folder, self.rootname + "_c0f"))
+            stokes = stokes.data
+        else:
+            raise ValueError("Input must be a path to a fits file or an HDUlist")
+        # FOS spectra are given in flux density with respect to angstrom wavelength
+        self.hd["DENSITY"] = True
+        self.hd["XUNIT"], self.hd["YUNIT"] = r"Wavelength [$\AA$]", r"{0:s}F$_\lambda$ [$10^{{{1:d}}}$ erg s$^{{-1}}$ cm$^{{-2}} \AA^{{-1}}$]"
+
+        # We set the error to be half the distance to the next mesure
+        self.wav = np.concat((wav[0:2, :], wav[0].reshape(1, wav.shape[1]), wav[0].reshape(1, wav.shape[1])), axis=0)
+        self.wav_err = np.zeros((self.wav.shape[0], self.wav.shape[1], 2))
+        for i in range(1, self.wav.shape[1] - 1):
+            self.wav_err[:, i] = np.abs(
+                np.array([((self.wav[j][i] - self.wav[j][i - 1]) / 2.0, (self.wav[j][i + 1] - self.wav[j][i - 1]) / 2.0) for j in range(self.wav.shape[0])])
+            )
+        self.wav_err[:, 0] = np.array([self.wav_err[:, 1, 0], self.wav_err[:, 1, 0]]).T
+        self.wav_err[:, -1] = np.array([self.wav_err[:, -2, 1], self.wav_err[:, -2, 1]]).T
+
+        self.hd["MINWAV"], self.hd["MAXWAV"] = self.wav.min(), self.wav.max()
+        self.hd["STEPWAV"] = np.mean(self.wav_err) * 2.0
+        self.bin_edges = np.array(
+            [np.concat((self.wav[i] - self.wav_err[i, :, 0], [self.wav[i, -1] + self.wav_err[i, -1, -1]]), axis=0) for i in range(self.wav.shape[0])]
+        )
+
+        self.IQUV_cov = np.zeros((4, 4, self.wav.shape[0], self.wav.shape[1]))
+
+        # Special way of reading FOS spectropolarimetry fits files
+        self.I = stokes[0::14]
+        self.IQUV_cov[0, 0] = stokes[4::14] ** 2
+        self.Q = stokes[1::14]
+        self.IQUV_cov[1, 1] = stokes[5::14] ** 2
+        self.U = stokes[2::14]
+        self.IQUV_cov[2, 2] = stokes[6::14] ** 2
+        self.V = stokes[3::14]
+        self.IQUV_cov[3, 3] = stokes[7::14] ** 2
+        self.hd["DATAMIN"], self.hd["DATAMAX"] = self.I.min(), self.I.max()
+
+        # Each file contain 4 spectra: Pass 1, Pass 2, combination of the 2, Combination corrected for orientation and background
+        self.subspec = {}
+        for i, name in enumerate(["PASS1", "PASS2", "PASS12", "PASS12corr"]):
+            spec = specpol(self.wav[i].shape[0])
+            spec.hd, spec.wav, spec.wav_err, spec.I, spec.Q, spec.U, spec.V = self.hd, self.wav[i], self.wav_err[i], self.I[i], self.Q[i], self.U[i], self.V[i]
+            spec.bin_edges = np.concat((spec.wav - spec.wav_err[:, 0], [spec.wav[-1] + spec.wav_err[-1, 1]]), axis=0)
+            spec.hd["MINWAV"], spec.hd["MAXWAV"] = spec.wav.min(), spec.wav.max()
+            spec.hd["DATAMIN"], spec.hd["DATAMAX"] = spec.I.min(), spec.I.max()
+            spec.IQUV_cov = self.IQUV_cov[:, :, i, :]
+            # Only PASS12corr is corrected for telescope orientation
+            spec.rotate(-(name[-4:] != "corr") * spec.hd["PA_APER"])
+            self.subspec[name] = spec
+
+        # Following lines contain the polarization components computed by calfos
+        self.P_fos = stokes[8::14]
+        self.P_fos_err = stokes[11::14]
+        self.PC_fos = stokes[9::14]
+        self.PC_fos_err = stokes[12::14]
+        self.PA_fos = princ_angle(
+            np.degrees(stokes[10::14]) + np.concat((np.ones((3, stokes.shape[1])), np.zeros((1, stokes.shape[1]))), axis=0) * self.hd["PA_APER"]
+        )
+        self.PA_fos_err = princ_angle(np.degrees(stokes[13::14]))
+
+    def dump_txt(self, filename, spec_list=None, output_dir=""):
+        """
+        Dump current spectra to a text file.
+        """
+        outfiles = []
+        if spec_list is None:
+            spec_list = self.subspec
+        for i, name in enumerate(["PASS1", "PASS2", "PASS12", "PASS12corr"]):
+            outfiles.append(spec_list[name].dump_txt(filename="_".join([filename, name]), output_dir=output_dir))
+        return outfiles
+
+    def plot(self, spec_list=None, savename=None, plots_folder="", fos=False):
+        """
+        Display current spectra in single figure.
+        """
+        outfiles = []
+        if hasattr(self, "ax"):
+            del self.ax
+        if hasattr(self, "fig"):
+            del self.fig
+        if spec_list is None:
+            spec_list = self.subspec
+        self.fig, self.ax = plt.subplots(4, 2, sharex=True, sharey="col", figsize=(20, 10), layout="constrained")
+        for i, (name, title) in enumerate(
+            [("PASS1", "Pass Direction 1"), ("PASS2", "Pass Direction 2"), ("PASS12", "Pass Direction 1&2"), ("PASS12corr", "Pass Direction 1&2 corrected")]
+        ):
+            self.ax[i][0].set_title(title)
+            if fos:
+                self.ax[i][0] = spec_list[name].plot(ax=self.ax[i][0])
+                self.ax[i][1].set_xlabel(r"Wavelength [$\AA$]")
+                secax1 = self.ax[i][0].secondary_xaxis("top", functions=(self.rest, self.unrest))
+                secax1.set_xlabel(r"Rest wavelength [$\AA$]")
+                secax2 = self.ax[i][1].secondary_xaxis("top", functions=(self.rest, self.unrest))
+                secax2.set_xlabel(r"Rest wavelength [$\AA$]")
+                self.ax[i][1].errorbar(self.wav[i], self.P_fos[i], xerr=self.wav_err[i].T, yerr=self.P_fos_err[i], color="b", fmt=".", label="P_fos")
+                self.ax[i][1].set_ylim([0.0, 1.0])
+                self.ax[i][1].set_ylabel(r"P", color="b")
+                self.ax[i][1].tick_params(axis="y", color="b", labelcolor="b")
+                ax22 = self.ax[i][1].twinx()
+                ax22.errorbar(self.wav[i], self.PA_fos[i], xerr=self.wav_err[i].T, yerr=self.PA_fos_err[i], color="r", fmt=".", label="PA_fos [°]")
+                ax22.set_ylim([0.0, 180.0])
+                ax22.set_ylabel(r"PA", color="r")
+                ax22.tick_params(axis="y", color="r", labelcolor="r")
+                h2, l2 = self.ax[i][1].get_legend_handles_labels()
+                h22, l22 = ax22.get_legend_handles_labels()
+                self.ax[i][1].legend(h2 + h22, l2 + l22, ncols=2, loc=1)
+            else:
+                self.ax[i] = spec_list[name].plot(ax=self.ax[i])
+        # self.ax[0][0].set_ylim(ymin=0.0)
+
+        self.fig.suptitle("_".join([self.hd["TARGNAME"], str(self.hd["PROPOSID"]), self.hd["ROOTNAME"], self.hd["APER_ID"]]))
+        if savename is not None:
+            self.fig.savefig(join_path(plots_folder, savename + ".pdf"), dpi=300, bbox_inches="tight")
+            outfiles.append(join_path(plots_folder, savename + ".pdf"))
+        return outfiles
+
+    def bin_size(self, size):
+        """
+        Rebin spectra to selected bin size in Angstrom.
+        """
+        key = "{0:.2f}bin".format(size)
+        if key not in self.subspec.keys():
+            self.subspec[key] = dict([])
+            for name in ["PASS1", "PASS2", "PASS12", "PASS12corr"]:
+                self.subspec[key][name] = self.subspec[name].bin_size(size)
+        return self.subspec[key]
+
+    def __add__(self, other):
+        """
+        Spectra addition, if not same binning default to self bins.
+        """
+        spec_a = FOSspecpol(self)
+        if np.all(self.wav == other.wav):
+            spec_b = other
+        else:
+            bin_edges = np.zeros(spec_a.wav.shape[0] + 1)
+            bin_edges[:-1], bin_edges[-1] = spec_a.wav - spec_a.wav_err[:, 0], spec_a.wav[-1] + spec_a.wav_err[-1:1]
+            spec_b = other.bin(bin_edges=bin_edges)
+
+        spec_a.I += deepcopy(spec_b.I)
+        spec_a.Q += deepcopy(spec_b.Q)
+        spec_a.U += deepcopy(spec_b.U)
+        spec_a.V += deepcopy(spec_b.V)
+        spec_a.IQUV_cov += deepcopy(spec_b.IQUV_cov)
+        for name in ["PASS1", "PASS2", "PASS12", "PASS12corr"]:
+            spec_a.subspec[name] += deepcopy(spec_b.subspec[name])
+            spec_a.subspec[name].hd["DATAMIN"], spec_a.subspec[name].hd["DATAMAX"] = spec_a.subspec[name].I.min(), spec_a.subspec[name].I.max()
+            spec_a.subspec[name].hd["EXPTIME"] += spec_b.subspec[name].hd["EXPTIME"]
+            spec_a.subspec[name].hd["ROOTNAME"] += "+" + spec_b.subspec[name].hd["ROOTNAME"]
+
+        spec_a.hd["DATAMIN"], spec_a.hd["DATAMAX"] = spec_a.I.min(), spec_a.I.max()
+        spec_a.hd["EXPTIME"] += spec_b.hd["EXPTIME"]
+        spec_a.hd["ROOTNAME"] += "+" + spec_b.hd["ROOTNAME"]
+        return spec_a
+
+    def __deepcopy__(self, memo):
+        spec = FOSspecpol(self.wav.shape[0])
+        spec.__dict__.update(self.__dict__)
+
+        for key in self.subspec.keys():
+            spec.subspec[key] = deepcopy(self.subspec[key])
+
+        return spec
+
+    def __del__(self):
+        if hasattr(self, "ax"):
+            del self.ax
+        if hasattr(self, "fig"):
+            del self.fig
+
+
+def main(infiles, bin_size=None, output_dir=None):
+    """
+    Produce (binned and summed) spectra for a list of given fits files.
+    """
+    outfiles = []
+    if infiles is not None:
+        # Divide path in folder + filename
+        prod = np.array([["/".join(filepath.split("/")[:-1]), filepath.split("/")[-1]] for filepath in infiles], dtype=str)
+        obs_dir = np.unique(["/".join(file.split("/")[:-1]) for file in infiles])
+        for dir in obs_dir:
+            # Create missing data/plot folder for tydiness
+            if not path_exists(dir):
+                system("mkdir -p {0:s} {1:s}".format(dir, dir.replace("data", "plots")))
+    else:
+        print("Must input files to process.")
+        return 1
+
+    data_folder = np.unique(prod[:, 0])
+    if output_dir is None:
+        output_dir = data_folder[0]
+    try:
+        plots_folder = output_dir.replace("data", "plots")
+    except ValueError:
+        plots_folder = output_dir
+    if not path_exists(plots_folder):
+        system("mkdir -p {0:s} ".format(plots_folder))
+
+    aper = dict([])
+    roots = np.unique([p[1].split("_")[0] for p in prod])
+    # Iteration on each observation in infiles
+    for rootname in roots:
+        print(rootname)
+        if data_folder.shape[0] > 1:
+            # For multiple folders (multiple filters) match data_folder on file rootname
+            spec = FOSspecpol(rootname, prod[np.array([p[1].split("_")[0] == rootname for p in prod])][0, 0])
+        else:
+            spec = FOSspecpol(rootname, data_folder[0])
+        filename = "_".join([spec.hd["TARGNAME"], "FOS", str(spec.hd["PROPOSID"]), spec.rootname, spec.hd["APER_ID"]])
+        if bin_size is not None:
+            key = "{0:.2f}bin".format(bin_size)
+            spec.bin_size(bin_size)
+            # Only output binned spectra
+            outfiles += spec.dump_txt("_".join([filename, key]), spec_list=spec.subspec[key], output_dir=output_dir)
+            outfiles += spec.plot(savename="_".join([filename, key]), spec_list=spec.subspec[key], plots_folder=plots_folder)
+
+            # Save corrected and combined pass for later summation, only sum on same aperture
+            if spec.hd["APER_ID"] in aper.keys():
+                aper[str(spec.hd["APER_ID"])].append(specpol(spec.subspec[key]["PASS12corr"]))
+            else:
+                aper[str(spec.hd["APER_ID"])] = [specpol(spec.subspec[key]["PASS12corr"])]
+        else:
+            outfiles += spec.dump_txt(filename, output_dir=output_dir)
+            outfiles += spec.plot(savename=filename, plots_folder=plots_folder)
+            if spec.hd["APER_ID"] in aper.keys():
+                aper[str(spec.hd["APER_ID"])].append(specpol(spec.subspec["PASS12corr"]))
+            else:
+                aper[str(spec.hd["APER_ID"])] = [specpol(spec.subspec["PASS12corr"])]
+    plt.close("all")
+
+    # Sum spectra acquired through same aperture
+    for key in aper.keys():
+        rootnames = [s.hd["ROOTNAME"] for s in aper[key]]
+        print(*rootnames)
+        spec = np.sum(aper[key])
+        spec.hd["ROOTNAME"] = "".join(p for p, *r in zip(*rootnames) if all(p == c for c in r)) + "_SUM"
+        filename = "_".join([spec.hd["TARGNAME"], "FOS", str(spec.hd["PROPOSID"]), spec.hd["ROOTNAME"]])
+        if bin_size is not None:
+            filename += "_{0:.2f}bin".format(bin_size)
+        # Output summed spectra
+        outfiles.append(spec.dump_txt("_".join([filename, key]), output_dir=output_dir))
+        outfiles.append(spec.plot(savename="_".join([filename, key]), plots_folder=plots_folder)[2])
+    plt.show()
+
+    return outfiles
+
+
+if __name__ == "__main__":
+    import argparse
+
+    parser = argparse.ArgumentParser(description="Display and dump FOS Spectropolarimetry")
+    parser.add_argument("-f", "--files", metavar="path", required=False, nargs="*", help="the full or relative path to the data products", default=None)
+    parser.add_argument("-b", "--bin", metavar="bin_size", required=False, help="The bin size to resample spectra", type=float, default=None)
+    parser.add_argument(
+        "-o", "--output_dir", metavar="directory_path", required=False, help="output directory path for the data products", type=str, default=None
+    )
+    args = parser.parse_args()
+    exitcode = main(infiles=args.files, bin_size=args.bin, output_dir=args.output_dir)
+    print("Written to: ", exitcode)

requirements.txt

@@ -0,0 +1,7 @@
+# Requirements for FOC_Reduction
+
+numpy
+scipy
+astropy
+astroquery  # To directly query the Mast archives
+matplotlib >= 3.8   # Brought some depreciation on ContourSet