diff --git a/iop4lib/db/photopolresult.py b/iop4lib/db/photopolresult.py index 2e16eeef..db83762b 100644 --- a/iop4lib/db/photopolresult.py +++ b/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 diff --git a/iop4lib/instruments/cafos.py b/iop4lib/instruments/cafos.py index e341815c..b00f7db4 100644 --- a/iop4lib/instruments/cafos.py +++ b/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 diff --git a/iop4lib/instruments/dipol.py b/iop4lib/instruments/dipol.py index d4ff6443..2d13f5d8 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, +) # 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 diff --git a/iop4lib/instruments/instrument.py b/iop4lib/instruments/instrument.py index e0dd2638..137789c2 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("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=}") diff --git a/iop4lib/utils/polarization.py b/iop4lib/utils/polarization.py index 48bf3180..472876ea 100644 --- a/iop4lib/utils/polarization.py +++ b/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']):