From 7dab4e38b3e7865f474086e96fcea1105ccf528e Mon Sep 17 00:00:00 2001 From: Juan Escudero Date: Sat, 4 Jul 2026 23:20:35 +0000 Subject: [PATCH 1/6] feat: implement 1-pair relative stokes fit --- iop4lib/instruments/cafos.py | 12 +++- iop4lib/instruments/dipol.py | 7 ++- iop4lib/instruments/instrument.py | 29 +++++----- iop4lib/utils/polarization.py | 91 +++++++++++++++++++++++++++++++ 4 files changed, 120 insertions(+), 19 deletions(-) diff --git a/iop4lib/instruments/cafos.py b/iop4lib/instruments/cafos.py index e341815c..e1449199 100644 --- a/iop4lib/instruments/cafos.py +++ b/iop4lib/instruments/cafos.py @@ -20,9 +20,14 @@ from .instrument import InstrumentHWP from iop4lib.telescopes import CAHAT220 from iop4lib.utils.polarization import ( + # compute_stokes_HWP_fit_full # compute_stokes_HWP_analytical, compute_stokes_HWP_fit_rel, # compute_stokes_HWP_fit_rel_nonideal, + # compute_stokes_HWP_fit_1pair + # compute_stokes_HWP_analytical_1pair, + compute_stokes_HWP_fit_1pair_rel, + # compute_stokes_HWP_fit_1pair_rel_nonideal, ) # logging @@ -61,9 +66,10 @@ 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_full, compute_stokes_HWP_fit_1pair + # default_pol_method = compute_stokes_HWP_analytical, compute_stokes_HWP_analytical_1pair + default_pol_method = compute_stokes_HWP_fit_rel, compute_stokes_HWP_fit_1pair_rel + # default_pol_method = compute_stokes_HWP_fit_rel_nonideal, compute_stokes_HWP_fit_1pair_rel_nonideal rot_angles_expected = {0.0, 22.48, 44.98, 67.48} min_rot_angles_required = 4 diff --git a/iop4lib/instruments/dipol.py b/iop4lib/instruments/dipol.py index d4ff6443..6cf6a779 100644 --- a/iop4lib/instruments/dipol.py +++ b/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_nonideal, +) # 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_nonideal 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 diff --git a/iop4lib/instruments/instrument.py b/iop4lib/instruments/instrument.py index e0dd2638..4db87fcf 100644 --- a/iop4lib/instruments/instrument.py +++ b/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("dF0") - 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=}") diff --git a/iop4lib/utils/polarization.py b/iop4lib/utils/polarization.py index 48bf3180..4639838c 100644 --- a/iop4lib/utils/polarization.py +++ b/iop4lib/utils/polarization.py @@ -1153,6 +1153,97 @@ 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) + dE = 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, + ) + +@polmethod(name="HWP_fit_1pair_rel_nonideal") +def compute_stokes_HWP_fit_1pair_rel_nonideal( + 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 (relative, non-ideal HWP).""" + + 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_nonideal( + theta=theta, + FO=FO, + FE=FE, + dFO=dFO, + dFE=dFE, + inst_pol_dict=inst_pol_dict, + plot=plot, + fig=fig, + annotate=annotate, + ) class PolarimetryGroup(list['ReducedFit']): From 4509011ea111df7cd4af8291c91f85c21e23b680 Mon Sep 17 00:00:00 2001 From: Juan Escudero Date: Sat, 4 Jul 2026 23:40:09 +0000 Subject: [PATCH 2/6] fix photopolresult plotting --- iop4lib/db/photopolresult.py | 77 ++++++++++++++++++++++-------------- 1 file changed, 47 insertions(+), 30 deletions(-) diff --git a/iop4lib/db/photopolresult.py b/iop4lib/db/photopolresult.py index 2e16eeef..2cb0d7ff 100644 --- a/iop4lib/db/photopolresult.py +++ b/iop4lib/db/photopolresult.py @@ -575,10 +575,14 @@ 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, + compute_stokes_HWP_fit_1pair_rel_nonideal, ) fig = fig or plt.gcf() @@ -607,19 +611,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("dF0") - 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 +644,29 @@ 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, + compute_stokes_HWP_fit_1pair_rel_nonideal, + ] - 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" @@ -679,8 +697,7 @@ 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 all_models = kwargs.get("all_models", False) default_width = default_width if all_models else 1024 From becb931e94d7dae7fe08a992bb3e4a57bb02426c Mon Sep 17 00:00:00 2001 From: Juan Escudero Date: Sun, 5 Jul 2026 17:55:45 +0000 Subject: [PATCH 3/6] fix --- iop4lib/db/photopolresult.py | 2 +- iop4lib/instruments/instrument.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/iop4lib/db/photopolresult.py b/iop4lib/db/photopolresult.py index 2cb0d7ff..406d4e2c 100644 --- a/iop4lib/db/photopolresult.py +++ b/iop4lib/db/photopolresult.py @@ -625,7 +625,7 @@ def plot_polarimetry(self, fig=None, all_models=False): args_dict.pop("dFE") elif only_pair == "E": args_dict.pop("FO") - args_dict.pop("dF0") + args_dict.pop("dFO") if not all_models: diff --git a/iop4lib/instruments/instrument.py b/iop4lib/instruments/instrument.py index 4db87fcf..137789c2 100644 --- a/iop4lib/instruments/instrument.py +++ b/iop4lib/instruments/instrument.py @@ -1321,7 +1321,7 @@ def compute_relative_polarimetry(cls, polarimetry_group: 'PolarimetryGroup') -> args_dict.pop("dFE") elif only_pair == "E": args_dict.pop("FO") - args_dict.pop("dF0") + args_dict.pop("dFO") 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})" From 71b9ec6a41ddcdb219c495301fc42227e0ebf265 Mon Sep 17 00:00:00 2001 From: Juan Escudero Date: Sun, 5 Jul 2026 20:26:42 +0000 Subject: [PATCH 4/6] fix --- iop4lib/utils/polarization.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/iop4lib/utils/polarization.py b/iop4lib/utils/polarization.py index 4639838c..26a72597 100644 --- a/iop4lib/utils/polarization.py +++ b/iop4lib/utils/polarization.py @@ -1197,7 +1197,7 @@ def compute_stokes_HWP_fit_1pair_rel( if FO is not None: FE = np.roll(FO, 1) - dE = np.roll(dFO, 1) + dFE = np.roll(dFO, 1) else: FO = np.roll(FE, 1) dFO = np.roll(dFE, 1) From a207fca9672c25727880f9f9cfd405f9f115a23d Mon Sep 17 00:00:00 2001 From: Juan Escudero Date: Sun, 5 Jul 2026 20:46:19 +0000 Subject: [PATCH 5/6] remove HWP_fit_1pair_rel_nonideal Kappa does not make sense (1pair uses one beam only). --- iop4lib/instruments/cafos.py | 9 --------- iop4lib/instruments/dipol.py | 4 ++-- iop4lib/utils/polarization.py | 28 ---------------------------- 3 files changed, 2 insertions(+), 39 deletions(-) diff --git a/iop4lib/instruments/cafos.py b/iop4lib/instruments/cafos.py index e1449199..b00f7db4 100644 --- a/iop4lib/instruments/cafos.py +++ b/iop4lib/instruments/cafos.py @@ -20,14 +20,8 @@ from .instrument import InstrumentHWP from iop4lib.telescopes import CAHAT220 from iop4lib.utils.polarization import ( - # compute_stokes_HWP_fit_full - # compute_stokes_HWP_analytical, compute_stokes_HWP_fit_rel, - # compute_stokes_HWP_fit_rel_nonideal, - # compute_stokes_HWP_fit_1pair - # compute_stokes_HWP_analytical_1pair, compute_stokes_HWP_fit_1pair_rel, - # compute_stokes_HWP_fit_1pair_rel_nonideal, ) # logging @@ -66,10 +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_fit_full, compute_stokes_HWP_fit_1pair - # default_pol_method = compute_stokes_HWP_analytical, compute_stokes_HWP_analytical_1pair default_pol_method = compute_stokes_HWP_fit_rel, compute_stokes_HWP_fit_1pair_rel - # default_pol_method = compute_stokes_HWP_fit_rel_nonideal, compute_stokes_HWP_fit_1pair_rel_nonideal rot_angles_expected = {0.0, 22.48, 44.98, 67.48} min_rot_angles_required = 4 diff --git a/iop4lib/instruments/dipol.py b/iop4lib/instruments/dipol.py index 6cf6a779..2d13f5d8 100644 --- a/iop4lib/instruments/dipol.py +++ b/iop4lib/instruments/dipol.py @@ -33,7 +33,7 @@ from iop4lib.telescopes import OSNT090 from iop4lib.utils.polarization import ( compute_stokes_HWP_fit_rel_nonideal, - compute_stokes_HWP_fit_1pair_rel_nonideal, + compute_stokes_HWP_fit_1pair_rel, ) # logging @@ -76,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, compute_stokes_HWP_fit_1pair_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 diff --git a/iop4lib/utils/polarization.py b/iop4lib/utils/polarization.py index 26a72597..472876ea 100644 --- a/iop4lib/utils/polarization.py +++ b/iop4lib/utils/polarization.py @@ -1215,35 +1215,7 @@ def compute_stokes_HWP_fit_1pair_rel( kappa=kappa, ) -@polmethod(name="HWP_fit_1pair_rel_nonideal") -def compute_stokes_HWP_fit_1pair_rel_nonideal( - 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 (relative, non-ideal HWP).""" - - 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_nonideal( - theta=theta, - FO=FO, - FE=FE, - dFO=dFO, - dFE=dFE, - inst_pol_dict=inst_pol_dict, - plot=plot, - fig=fig, - annotate=annotate, - ) class PolarimetryGroup(list['ReducedFit']): From deb4787978763f5abfd3a425ec502c3f5b279908 Mon Sep 17 00:00:00 2001 From: Juan Escudero Date: Sun, 5 Jul 2026 20:53:48 +0000 Subject: [PATCH 6/6] fix --- iop4lib/db/photopolresult.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/iop4lib/db/photopolresult.py b/iop4lib/db/photopolresult.py index 406d4e2c..db83762b 100644 --- a/iop4lib/db/photopolresult.py +++ b/iop4lib/db/photopolresult.py @@ -582,7 +582,6 @@ def plot_polarimetry(self, fig=None, all_models=False): compute_stokes_HWP_fit_1pair, compute_stokes_HWP_analytical_1pair, compute_stokes_HWP_fit_1pair_rel, - compute_stokes_HWP_fit_1pair_rel_nonideal, ) fig = fig or plt.gcf() @@ -650,7 +649,6 @@ def plot_polarimetry(self, fig=None, all_models=False): compute_stokes_HWP_fit_1pair, compute_stokes_HWP_fit_1pair_rel, compute_stokes_HWP_analytical_1pair, - compute_stokes_HWP_fit_1pair_rel_nonideal, ] gs = fig.add_gridspec(1, len(pol_methods)) @@ -697,8 +695,9 @@ def plot_polarimetry(self, fig=None, all_models=False): def get_img_polarimetry(self, **kwargs): - default_width = 4*1024 - + only_pair = self.astrosource.metadata.get(f"{self.instrument}.polarimetry.only_pair") + 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