feat: implement 1-pair relative stokes fit

iop4lib/db/photopolresult.py

@@ -575,10 +575,13 @@ def plot_polarimetry(self, fig=None, all_models=False):
         from iop4lib.utils import get_column_values
         from iop4lib.utils.polarization import (
             compute_stokes_HWP_fit_full,
-            compute_stokes_HWP_fit_rel,
-            compute_stokes_HWP_fit_1pair,
             compute_stokes_HWP_analytical,
+            compute_stokes_HWP_fit_rel,
             compute_stokes_HWP_fit_rel_nonideal,
+            # and equivalent for 1 pair only
+            compute_stokes_HWP_fit_1pair,
+            compute_stokes_HWP_analytical_1pair,
+            compute_stokes_HWP_fit_1pair_rel,
         )
 
         fig = fig or plt.gcf()
@@ -607,19 +610,31 @@ def plot_polarimetry(self, fig=None, all_models=False):
         redf_0 = self.reducedfits.first()
         inst_pol_dict = redf_0.instrument_cls.get_instrumental_polarization(redf_0)
 
-        if not (only_pair := astrosource.metadata.get(f"{self.instrument}.polarimetry.only_pair")):
+        only_pair = astrosource.metadata.get(f"{self.instrument}.polarimetry.only_pair")
 
-            if not all_models:
+        args_dict = dict(
+            theta=theta,
+            FO=FO, dFO=dFO,
+            FE=FE, dFE=dFE,
+            inst_pol_dict=inst_pol_dict,
+        )
 
-                gs = fig.add_gridspec(1, 1)
-                subfig = fig.add_subfigure(gs[0,0])
+        if only_pair == "O":
+            args_dict.pop("FE")
+            args_dict.pop("dFE")
+        elif only_pair == "E":
+            args_dict.pop("FO")
+            args_dict.pop("dFO")
 
-                pol_method = self.instrument_cls.default_pol_method
-                stokes_nocorr, fit_stats = pol_method(theta, FO=FO, dFO=dFO, FE=FE, dFE=dFE, inst_pol_dict=inst_pol_dict, plot=True, annotate=True, fig=subfig)
-                title = subfig.suptitle(pol_method.name, y=1)
-                method_name = pol_method.name
+        if not all_models:
 
-            else:
+            pol_methods = [
+                self.instrument_cls.default_pol_method[0 if not only_pair else 1]
+            ]
+
+        else:
+
+            if not only_pair:
 
                 pol_methods = [
                     compute_stokes_HWP_fit_full,
@@ -628,27 +643,28 @@ def plot_polarimetry(self, fig=None, all_models=False):
                     compute_stokes_HWP_fit_rel_nonideal,
                 ]
 
-                gs = fig.add_gridspec(1, len(pol_methods))
+            else:
 
-                for i, pol_method in enumerate(pol_methods):
+                pol_methods = [
+                    compute_stokes_HWP_fit_1pair,
+                    compute_stokes_HWP_fit_1pair_rel,
+                    compute_stokes_HWP_analytical_1pair,
+                ]
 
-                    subfig = fig.add_subfigure(gs[0,i])
+        gs = fig.add_gridspec(1, len(pol_methods))
 
-                    stokes_nocorr_i, fit_stats_i = pol_method(theta, FO=FO, dFO=dFO, FE=FE, dFE=dFE, inst_pol_dict=inst_pol_dict, plot=True, annotate=True, fig=subfig)
-                    title = subfig.suptitle(pol_method.name)
+        for i, pol_method in enumerate(pol_methods):
 
-                    if pol_method == self.instrument_cls.default_pol_method:
-                        title.set(color='green', weight='bold')
-                        method_name = pol_method.name
-                        stokes_nocorr = stokes_nocorr_i
-
-        else:
-            if only_pair == 'O':
-                method_name = "compute_stokes_HWP_fit_1pair (O)"
-                stokes_nocorr, fit_stats = compute_stokes_HWP_fit_1pair(theta, FO=FO, dFO=dFO, inst_pol_dict=inst_pol_dict, plot=True, annotate=True, fig=fig)
-            elif only_pair == 'E':
-                method_name = "compute_stokes_HWP_fit_1pair (E)"
-                stokes_nocorr, fit_stats = compute_stokes_HWP_fit_1pair(theta, FE=FE, dFE=dFE, inst_pol_dict=inst_pol_dict, plot=True, annotate=True, fig=fig)
+            subfig = fig.add_subfigure(gs[0,i])
+
+            args_dict_i = args_dict | dict(plot=True, annotate=True, fig=subfig)
+            stokes_nocorr_i, fit_stats_i = pol_method(**args_dict_i)
+            title = subfig.suptitle(pol_method.name)
+
+            if pol_method == self.instrument_cls.default_pol_method[0 if not only_pair else 1]:
+                title.set(color='green', weight='bold')
+                method_name =  pol_method.name if not only_pair else f"{pol_method.name} ({only_pair})"
+                stokes_nocorr = stokes_nocorr_i
 
         fig_title = (
             f"{self}\n"
@@ -680,8 +696,8 @@ def plot_polarimetry(self, fig=None, all_models=False):
     def get_img_polarimetry(self, **kwargs):
 
         only_pair = self.astrosource.metadata.get(f"{self.instrument}.polarimetry.only_pair")
-        default_width = 4*1024 if not only_pair else 1024
-
+        default_width = 4*1024 if not only_pair else 3*1024
+        
         all_models = kwargs.get("all_models", False)
         default_width = default_width if all_models else 1024
 

iop4lib/instruments/cafos.py

@@ -20,9 +20,8 @@
 from .instrument import InstrumentHWP
 from iop4lib.telescopes import CAHAT220
 from iop4lib.utils.polarization import (
-    # compute_stokes_HWP_analytical,
     compute_stokes_HWP_fit_rel,
-    # compute_stokes_HWP_fit_rel_nonideal,
+    compute_stokes_HWP_fit_1pair_rel,
 )
 
 # logging
@@ -61,9 +60,7 @@ class CAFOS(InstrumentHWP):
     # TODO: which polarimetry method is better will depend on the instrument
     # and the number of observations.
 
-    # default_pol_method = compute_stokes_HWP_analytical
-    default_pol_method = compute_stokes_HWP_fit_rel
-    # default_pol_method = compute_stokes_HWP_fit_rel_nonideal
+    default_pol_method = compute_stokes_HWP_fit_rel, compute_stokes_HWP_fit_1pair_rel
 
     rot_angles_expected = {0.0, 22.48, 44.98, 67.48}
     min_rot_angles_required = 4

iop4lib/instruments/dipol.py

@@ -31,7 +31,10 @@
 from iop4lib.utils.plotting import imshow_w_sources, plot_preview_astrometry
 from iop4lib.utils.astrometry import BuildWCSResult
 from iop4lib.telescopes import OSNT090
-from iop4lib.utils.polarization import compute_stokes_HWP_fit_rel_nonideal
+from iop4lib.utils.polarization import (
+    compute_stokes_HWP_fit_rel_nonideal,
+    compute_stokes_HWP_fit_1pair_rel,
+)
 
 # logging
 import logging
@@ -73,7 +76,7 @@ class DIPOL(InstrumentHWP):
     # In [21]: np.mean(angle_L), np.median(angle_L), np.std(angle_L)
     # Out[21]: (179.08971366048235, 177.6282921156412, 1.9341471640122656)
 
-    default_pol_method = compute_stokes_HWP_fit_rel_nonideal
+    default_pol_method = compute_stokes_HWP_fit_rel_nonideal, compute_stokes_HWP_fit_1pair_rel
 
     rot_angles_expected = {0.0, 22.5, 45.0, 67.5, 90.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0, 247.5, 270.0, 292.5, 315.0, 337.5}
     min_rot_angles_required = 8

iop4lib/instruments/instrument.py

@@ -46,9 +46,6 @@
     calibrate_photopolresult,
     NoCalibratorsFound,
 )
-from iop4lib.utils.polarization import (
-    compute_stokes_HWP_fit_1pair,
-)
 
 # logging
 import logging
@@ -1313,20 +1310,24 @@ def compute_relative_polarimetry(cls, polarimetry_group: 'PolarimetryGroup') ->
                 # read possible special case from target source metadata
                 only_pair = astrosource.metadata.get(f"{cls.name}.polarimetry.only_pair")
 
-                if not only_pair:
+                args_dict = dict(
+                    theta=theta,
+                    FO=FO, dFO=dFO,
+                    FE=FE, dFE=dFE,
+                )
 
-                    logger.info(f"Computing stokes parameters with {cls.default_pol_method.name}")
+                if only_pair == "O":
+                    args_dict.pop("FE")
+                    args_dict.pop("dFE")
+                elif only_pair == "E":
+                    args_dict.pop("FO")
+                    args_dict.pop("dFO")
 
-                    stokes_nocorr, fit_stats = cls.default_pol_method(theta, FO=FO, dFO=dFO, FE=FE, dFE=dFE)
-
-                else:
-
-                    logger.info(f"This source metadata indicates to use only the {only_pair} pair, computing stoke parameters with compute_stokes_HWP_fit_1pair")
+                pol_method = cls.default_pol_method[0 if not only_pair else 1]
+                method_name = pol_method.name if not only_pair else f"{pol_method.name} ({only_pair})"
 
-                    if only_pair == 'O':
-                        stokes_nocorr, fit_stats = compute_stokes_HWP_fit_1pair(theta, FO=FO, dFO=dFO)
-                    elif only_pair == 'E':
-                        stokes_nocorr, fit_stats = compute_stokes_HWP_fit_1pair(theta, FE=FE, dFE=dFE)
+                logger.info(f"Computing stokes parameters with {method_name}")
+                stokes_nocorr, fit_stats = pol_method(**args_dict)
 
                 logger.debug(f"{stokes_nocorr=}")
                 logger.debug(f"{fit_stats=}")

iop4lib/utils/polarization.py

@@ -1153,6 +1153,69 @@ def compute_stokes_HWP_fit_1pair(
     return stokes, fit_stats
 
 
+@polmethod(name="HWP_analytical_1pair")
+def compute_stokes_HWP_analytical_1pair(
+        theta, FO=None, dFO=None, FE=None, dFE=None,
+        inst_pol_dict=None,
+        plot=False, fig=None, annotate=False,
+    ):
+    """Compute polarimetry fitting only the O (or E) pair (analytical)."""
+
+    assert (FO is None) ^ bool(FE is None), "must specify one and only one of FO or FE"
+    assert bool(dFO is None) ^ bool(dFE is None), "must specify one and only one of dFO or dFE"
+
+    if FO is not None:
+        FE = np.roll(FO, 1)
+        dFE = np.roll(dFO, 1)
+    else:
+        FO = np.roll(FE, 1)
+        dFO = np.roll(dFE, 1)
+
+    return compute_stokes_HWP_analytical(
+        theta=theta,
+        FO=FO,
+        FE=FE,
+        dFO=dFO,
+        dFE=dFE,
+        inst_pol_dict=inst_pol_dict,
+        plot=plot,
+        fig=fig,
+        annotate=annotate,
+    )
+
+@polmethod(name="HWP_fit_1pair_rel")
+def compute_stokes_HWP_fit_1pair_rel(
+        theta, FO=None, dFO=None, FE=None, dFE=None,
+        inst_pol_dict=None,
+        plot=False, fig=None, annotate=False,
+        kappa=None,
+    ):
+    """Compute polarimetry fitting only the O (or E) pair (relative)."""
+
+    assert (FO is None) ^ bool(FE is None), "must specify one and only one of FO or FE"
+    assert bool(dFO is None) ^ bool(dFE is None), "must specify one and only one of dFO or dFE"
+
+    if FO is not None:
+        FE = np.roll(FO, 1)
+        dFE = np.roll(dFO, 1)
+    else:
+        FO = np.roll(FE, 1)
+        dFO = np.roll(dFE, 1)
+
+    return compute_stokes_HWP_fit_rel(
+        theta=theta,
+        FO=FO,
+        FE=FE,
+        dFO=dFO,
+        dFE=dFE,
+        inst_pol_dict=inst_pol_dict,
+        plot=plot,
+        fig=fig,
+        annotate=annotate,
+        kappa=kappa,
+    )
+
+
 
 class PolarimetryGroup(list['ReducedFit']):