From f4effac343d2059992adc16dfb8171bea7e8b784 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 8 Aug 2025 11:45:11 +0200 Subject: [PATCH] Add Stokes_cov_stat to fits and compute again P_debiased in plots --- package/FOC_reduction.py | 26 +++- package/lib/fits.py | 27 +++-- package/lib/plots.py | 254 ++++++++++++++++++++++++--------------- package/lib/reduction.py | 58 ++++----- 4 files changed, 223 insertions(+), 142 deletions(-) diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index 816283e..87dcbfa 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -41,7 +41,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Background estimation error_sub_type = "scott" # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) - subtract_error = 0.50 + subtract_error = 0.80 display_bkg = False # Data binning @@ -203,24 +203,38 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # FWHM of FOC have been estimated at about 0.03" across 1500-5000 Angstrom band, which is about 2 detector pixels wide # see Jedrzejewski, R.; Nota, A.; Hack, W. J., A Comparison Between FOC and WFPC2 # Bibcode : 1995chst.conf...10J - Stokes, Stokes_cov, header_stokes, s_IQU_stat = proj_red.compute_Stokes( + Stokes, Stokes_cov, header_stokes, Stokes_cov_stat = proj_red.compute_Stokes( data_array, error_array, data_mask, headers, FWHM=smoothing_FWHM, scale=smoothing_scale, smoothing=smoothing_function, transmitcorr=transmitcorr ) # Step 3: # Rotate images to have North up if rotate_North: - Stokes, Stokes_cov, data_mask, header_stokes, s_IQU_stat = proj_red.rotate_Stokes( - Stokes, Stokes_cov, data_mask, header_stokes, s_IQU_stat=s_IQU_stat, SNRi_cut=None + Stokes, Stokes_cov, data_mask, header_stokes, Stokes_cov_stat = proj_red.rotate_Stokes( + Stokes, Stokes_cov, data_mask, header_stokes, Stokes_cov_stat=Stokes_cov_stat, SNRi_cut=None ) # 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(Stokes, Stokes_cov, header_stokes, s_IQU_stat=s_IQU_stat) + P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(Stokes, Stokes_cov, header_stokes, Stokes_cov_stat=Stokes_cov_stat) # Step 4: # Save image to FITS. figname = "_".join([figname, figtype]) if figtype != "" else figname Stokes_hdul = proj_fits.save_Stokes( - 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 + Stokes, + Stokes_cov, + Stokes_cov_stat, + 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, ) outfiles.append("/".join([data_folder, Stokes_hdul[0].header["FILENAME"] + ".fits"])) diff --git a/package/lib/fits.py b/package/lib/fits.py index 1276b3f..b75a266 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -105,7 +105,9 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): return data_array, headers -def save_Stokes(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): +def save_Stokes( + Stokes, Stokes_cov, Stokes_cov_stat, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, header_stokes, data_mask, filename, data_folder="", return_hdul=False +): """ Save computed polarimetry parameters to a single fits file, updating header accordingly. @@ -116,8 +118,9 @@ def save_Stokes(Stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, Stokes parameters I, Q, U, V, Polarization degree and debieased, its error propagated and assuming Poisson noise, Polarization angle, its error propagated and assuming Poisson noise. - Stokes_cov : numpy.ndarray - Covariance matrix of the Stokes parameters I, Q, U. + Stokes_cov, Stokes_cov_stat : numpy.ndarray + Covariance matrix of the Stokes parameters I, Q, U, V and its statistical + uncertainties part. headers : header list Header of reference some keywords will be copied from (CRVAL, CDELT, INSTRUME, PROPOSID, TARGNAME, ORIENTAT, EXPTOT). @@ -135,14 +138,14 @@ def save_Stokes(Stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, ---------- Return: hdul : astropy.io.fits.hdu.hdulist.HDUList - HDUList containing the Stokes cube in the PrimaryHDU, then - P, s_P, PA, s_PA in this order. Headers have been updated to relevant - informations (WCS, orientation, data_type). + HDUList containing the Stokes cube in the PrimaryHDU, then Stokes_cov, + Stokes_cov_stat, P, s_P, PA, s_PA in this order. Headers have been updated + to relevant informations (WCS, orientation, data_type). Only returned if return_hdul is True. """ # Create new WCS object given the modified images new_wcs = WCS(header_stokes).celestial.deepcopy() - header = remove_stokes_axis_from_header(header_stokes).copy() + header = remove_stokes_axis_from_header(header_stokes).copy() if header_stokes["NAXIS"] > 2 else header_stokes.copy() if data_mask.shape != (1, 1): vertex = clean_ROI(data_mask) @@ -176,6 +179,7 @@ def save_Stokes(Stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, if data_mask.shape != (1, 1): Stokes = Stokes[:, vertex[2] : vertex[3], vertex[0] : vertex[1]] Stokes_cov = Stokes_cov[:, :, vertex[2] : vertex[3], vertex[0] : vertex[1]] + Stokes_cov_stat = Stokes_cov_stat[:, :, vertex[2] : vertex[3], vertex[0] : vertex[1]] P = P[vertex[2] : vertex[3], vertex[0] : vertex[1]] debiased_P = debiased_P[vertex[2] : vertex[3], vertex[0] : vertex[1]] s_P = s_P[vertex[2] : vertex[3], vertex[0] : vertex[1]] @@ -192,7 +196,7 @@ def save_Stokes(Stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, # Add I_stokes as PrimaryHDU header["datatype"] = ("STOKES", "type of data stored in the HDU") Stokes[np.broadcast_to((1 - data_mask).astype(bool), Stokes.shape)] = 0.0 - hdu_head = add_stokes_axis_to_header(header, 2) + hdu_head = add_stokes_axis_to_header(header, 0) primary_hdu = fits.PrimaryHDU(data=Stokes, header=hdu_head) primary_hdu.name = "STOKES" hdul.append(primary_hdu) @@ -200,6 +204,7 @@ def save_Stokes(Stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, # Add Stokes_cov, P, s_P, PA, s_PA to the HDUList for data, name in [ [Stokes_cov, "STOKES_COV"], + [Stokes_cov_stat, "STOKES_COV_STAT"], [P, "Pol_deg"], [debiased_P, "Pol_deg_debiased"], [s_P, "Pol_deg_err"], @@ -211,9 +216,9 @@ def save_Stokes(Stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, ]: hdu_head = header.copy() hdu_head["datatype"] = name - if name == "STOKES_COV": - hdu_head = add_stokes_axis_to_header(hdu_head, 2) - hdu_head = add_stokes_axis_to_header(hdu_head, 3) + if name[:10] == "STOKES_COV": + hdu_head = add_stokes_axis_to_header(hdu_head, 0) + hdu_head = add_stokes_axis_to_header(hdu_head, 0) data[np.broadcast_to((1 - data_mask).astype(bool), data.shape)] = 0.0 else: data[(1 - data_mask).astype(bool)] = 0.0 diff --git a/package/lib/plots.py b/package/lib/plots.py index 17ef000..2cfcdee 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -2113,7 +2113,7 @@ class crop_Stokes(crop_map): for dataset in self.hdul_crop: if dataset.header["datatype"] == "STOKES": dataset.data = deepcopy(dataset.data[:, vertex[2] : vertex[3], vertex[0] : vertex[1]]) - elif dataset.header["datatype"] == "STOKES_COV": + elif dataset.header["datatype"][:10] == "STOKES_COV": dataset.data = deepcopy(dataset.data[:, :, vertex[2] : vertex[3], vertex[0] : vertex[1]]) else: dataset.data = deepcopy(dataset.data[vertex[2] : vertex[3], vertex[0] : vertex[1]]) @@ -2142,6 +2142,12 @@ class crop_Stokes(crop_map): IQ_diluted_err = np.sqrt(np.sum(self.hdul_crop["STOKES_COV"].data[0, 1][mask] ** 2)) IU_diluted_err = np.sqrt(np.sum(self.hdul_crop["STOKES_COV"].data[0, 2][mask] ** 2)) QU_diluted_err = np.sqrt(np.sum(self.hdul_crop["STOKES_COV"].data[1, 2][mask] ** 2)) + I_diluted_stat_err = np.sqrt(np.sum(self.hdul_crop["STOKES_COV_STAT"].data[0, 0][mask])) + Q_diluted_stat_err = np.sqrt(np.sum(self.hdul_crop["STOKES_COV_STAT"].data[1, 1][mask])) + U_diluted_stat_err = np.sqrt(np.sum(self.hdul_crop["STOKES_COV_STAT"].data[2, 2][mask])) + IQ_diluted_stat_err = np.sqrt(np.sum(self.hdul_crop["STOKES_COV_STAT"].data[0, 1][mask] ** 2)) + IU_diluted_stat_err = np.sqrt(np.sum(self.hdul_crop["STOKES_COV_STAT"].data[0, 2][mask] ** 2)) + QU_diluted_stat_err = np.sqrt(np.sum(self.hdul_crop["STOKES_COV_STAT"].data[1, 2][mask] ** 2)) P_diluted = np.sqrt(Q_diluted**2 + U_diluted**2) / I_diluted P_diluted_err = (1.0 / I_diluted) * np.sqrt( @@ -2150,6 +2156,18 @@ class crop_Stokes(crop_map): - 2.0 * (Q_diluted / I_diluted) * IQ_diluted_err - 2.0 * (U_diluted / I_diluted) * IU_diluted_err ) + P_diluted_stat_err = ( + P_diluted + / I_diluted + * np.sqrt( + I_diluted_stat_err + - 2.0 / (I_diluted * P_diluted**2) * (Q_diluted * IQ_diluted_stat_err + U_diluted * IU_diluted_stat_err) + + 1.0 + / (I_diluted**2 * P_diluted**4) + * (Q_diluted**2 * Q_diluted_stat_err + U_diluted**2 * U_diluted_stat_err + 2.0 * Q_diluted * U_diluted * QU_diluted_stat_err) + ) + ) + debiased_P_diluted = np.sqrt(P_diluted**2 - P_diluted_stat_err**2) if P_diluted**2 > P_diluted_stat_err**2 else 0.0 PA_diluted = princ_angle((90.0 / np.pi) * np.arctan2(U_diluted, Q_diluted)) PA_diluted_err = (90.0 / (np.pi * (Q_diluted**2 + U_diluted**2))) * np.sqrt( @@ -2159,7 +2177,7 @@ class crop_Stokes(crop_map): for dataset in self.hdul_crop: if dataset.header["FILENAME"][-4:] != "crop": dataset.header["FILENAME"] += "_crop" - dataset.header["P_int"] = (P_diluted, "Integrated polarization degree") + dataset.header["P_int"] = (debiased_P_diluted, "Integrated polarization degree") dataset.header["sP_int"] = (np.ceil(P_diluted_err * 1000.0) / 1000.0, "Integrated polarization degree error") dataset.header["PA_int"] = (PA_diluted, "Integrated polarization angle") dataset.header["sPA_int"] = (np.ceil(PA_diluted_err * 10.0) / 10.0, "Integrated polarization angle error") @@ -2492,7 +2510,7 @@ class pol_map(object): ax_vec_sc = self.fig.add_axes([0.260, 0.030, 0.090, 0.01]) ax_snr_reset = self.fig.add_axes([0.060, 0.020, 0.05, 0.02]) ax_snr_conf = self.fig.add_axes([0.115, 0.020, 0.05, 0.02]) - SNRi_max = np.max(self.I[self.IQU_cov[0, 0] > 0.0] / np.sqrt(self.IQU_cov[0, 0][self.IQU_cov[0, 0] > 0.0])) + SNRi_max = np.max(self.I[self.Stokes_cov[0, 0] > 0.0] / np.sqrt(self.Stokes_cov[0, 0][self.Stokes_cov[0, 0] > 0.0])) SNRp_max = np.max(self.P[self.P_ERR > 0.0] / self.P_ERR[self.P_ERR > 0.0]) s_I_cut = Slider(ax_I_cut, r"$SNR^{I}_{cut}$", 1.0, int(SNRi_max * 0.95), valstep=0.5, valinit=self.SNRi_cut) self.P_ERR_cut = Slider( @@ -2988,9 +3006,13 @@ class pol_map(object): return self.U_ERR / np.where(self.I > 0, self.I, np.nan) @property - def IQU_cov(self): + def Stokes_cov(self): return self.Stokes["STOKES_COV"].data + @property + def Stokes_cov_stat(self): + return self.Stokes["STOKES_COV_STAT"].data + @property def P(self): return self.Stokes["POL_DEG_DEBIASED"].data @@ -3016,7 +3038,7 @@ class pol_map(object): @property def cut(self): - s_I = np.sqrt(self.IQU_cov[0, 0]) + s_I = np.sqrt(self.Stokes_cov[0, 0]) SNRp_mask, SNRi_mask = (np.zeros(self.P.shape).astype(bool), np.zeros(self.I.shape).astype(bool)) SNRi_mask[s_I > 0.0] = self.I[s_I > 0.0] / s_I[s_I > 0.0] > self.SNRi if self.SNRp >= 1.0: @@ -3159,7 +3181,7 @@ class pol_map(object): kwargs["alpha"] = 1.0 - 0.75 * (self.P < self.P_ERR) label = r"$\theta_{P}$ [°]" elif self.display_selection.lower() in ["snri"]: - s_I = np.sqrt(self.IQU_cov[0, 0]) + s_I = np.sqrt(self.Stokes_cov[0, 0]) SNRi = np.zeros(self.I.shape) SNRi[s_I > 0.0] = self.I[s_I > 0.0] / s_I[s_I > 0.0] self.data = SNRi @@ -3354,68 +3376,77 @@ class pol_map(object): ) fig.canvas.draw_idle() - def pol_int(self, fig=None, ax=None): + def pol_int(self, fig=None, ax=None, cut=False): str_conf = "" if self.region is None: - s_I = np.sqrt(self.IQU_cov[0, 0]) + s_I = np.sqrt(self.Stokes_cov[0, 0]) I_reg = self.I.sum() I_reg_err = np.sqrt(np.sum(s_I**2)) - P_reg = self.Stokes[0].header["P_int"] + debiased_P_reg = self.Stokes[0].header["P_int"] P_reg_err = self.Stokes[0].header["sP_int"] PA_reg = self.Stokes[0].header["PA_int"] PA_reg_err = self.Stokes[0].header["sPA_int"] - s_I = np.sqrt(self.IQU_cov[0, 0]) - s_Q = np.sqrt(self.IQU_cov[1, 1]) - s_U = np.sqrt(self.IQU_cov[2, 2]) - s_IQ = self.IQU_cov[0, 1] - s_IU = self.IQU_cov[0, 2] - s_QU = self.IQU_cov[1, 2] + if cut: + I_cut = self.I[self.cut].sum() + Q_cut = self.Q[self.cut].sum() + U_cut = self.U[self.cut].sum() + I_cut_err = np.sqrt(np.sum(self.Stokes_cov[0, 0][self.cut])) + Q_cut_err = np.sqrt(np.sum(self.Stokes_cov[1, 1][self.cut])) + U_cut_err = np.sqrt(np.sum(self.Stokes_cov[2, 2][self.cut])) + IQ_cut_err = np.sqrt(np.sum(self.Stokes_cov[0, 1][self.cut] ** 2)) + IU_cut_err = np.sqrt(np.sum(self.Stokes_cov[0, 2][self.cut] ** 2)) + QU_cut_err = np.sqrt(np.sum(self.Stokes_cov[1, 2][self.cut] ** 2)) + I_cut_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[0, 0][self.cut])) + Q_cut_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[1, 1][self.cut])) + U_cut_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[2, 2][self.cut])) + IQ_cut_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[0, 1][self.cut] ** 2)) + IU_cut_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[0, 2][self.cut] ** 2)) + QU_cut_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[1, 2][self.cut] ** 2)) - I_cut = self.I[self.cut].sum() - Q_cut = self.Q[self.cut].sum() - U_cut = self.U[self.cut].sum() - I_cut_err = np.sqrt(np.sum(s_I[self.cut] ** 2)) - Q_cut_err = np.sqrt(np.sum(s_Q[self.cut] ** 2)) - U_cut_err = np.sqrt(np.sum(s_U[self.cut] ** 2)) - IQ_cut_err = np.sqrt(np.sum(s_IQ[self.cut] ** 2)) - IU_cut_err = np.sqrt(np.sum(s_IU[self.cut] ** 2)) - QU_cut_err = np.sqrt(np.sum(s_QU[self.cut] ** 2)) - - with np.errstate(divide="ignore", invalid="ignore"): - P_cut = np.sqrt(Q_cut**2 + U_cut**2) / I_cut - P_cut_err = ( - np.sqrt( - (Q_cut**2 * Q_cut_err**2 + U_cut**2 * U_cut_err**2 + 2.0 * Q_cut * U_cut * QU_cut_err) / (Q_cut**2 + U_cut**2) - + ((Q_cut / I_cut) ** 2 + (U_cut / I_cut) ** 2) * I_cut_err**2 - - 2.0 * (Q_cut / I_cut) * IQ_cut_err - - 2.0 * (U_cut / I_cut) * IU_cut_err + with np.errstate(divide="ignore", invalid="ignore"): + P_cut = np.sqrt(Q_cut**2 + U_cut**2) / I_cut + P_cut_err = ( + np.sqrt( + (Q_cut**2 * Q_cut_err**2 + U_cut**2 * U_cut_err**2 + 2.0 * Q_cut * U_cut * QU_cut_err) / (Q_cut**2 + U_cut**2) + + ((Q_cut / I_cut) ** 2 + (U_cut / I_cut) ** 2) * I_cut_err**2 + - 2.0 * (Q_cut / I_cut) * IQ_cut_err + - 2.0 * (U_cut / I_cut) * IU_cut_err + ) + / I_cut ) - / I_cut - ) + P_cut_stat_err = ( + P_cut + / I_cut + * np.sqrt( + I_cut_stat_err + - 2.0 / (I_cut * P_cut**2) * (Q_cut * IQ_cut_stat_err + U_cut * IU_cut_stat_err) + + 1.0 / (I_cut**2 * P_cut**4) * (Q_cut**2 * Q_cut_stat_err + U_cut**2 * U_cut_stat_err + 2.0 * Q_cut * U_cut * QU_cut_stat_err) + ) + ) + debiased_P_cut = np.sqrt(P_cut**2 - P_cut_stat_err**2) if P_cut**2 > P_cut_stat_err**2 else 0.0 - PA_cut = princ_angle((90.0 / np.pi) * np.arctan2(U_cut, Q_cut)) - PA_cut_err = (90.0 / (np.pi * (Q_cut**2 + U_cut**2))) * np.sqrt( - U_cut**2 * Q_cut_err**2 + Q_cut**2 * U_cut_err**2 - 2.0 * Q_cut * U_cut * QU_cut_err - ) + PA_cut = princ_angle((90.0 / np.pi) * np.arctan2(U_cut, Q_cut)) + PA_cut_err = (90.0 / (np.pi * (Q_cut**2 + U_cut**2))) * np.sqrt( + U_cut**2 * Q_cut_err**2 + Q_cut**2 * U_cut_err**2 - 2.0 * Q_cut * U_cut * QU_cut_err + ) else: - s_I = np.sqrt(self.IQU_cov[0, 0]) - s_Q = np.sqrt(self.IQU_cov[1, 1]) - s_U = np.sqrt(self.IQU_cov[2, 2]) - s_IQ = self.IQU_cov[0, 1] - s_IU = self.IQU_cov[0, 2] - s_QU = self.IQU_cov[1, 2] - I_reg = self.I[self.region].sum() Q_reg = self.Q[self.region].sum() U_reg = self.U[self.region].sum() - I_reg_err = np.sqrt(np.sum(s_I[self.region] ** 2)) - Q_reg_err = np.sqrt(np.sum(s_Q[self.region] ** 2)) - U_reg_err = np.sqrt(np.sum(s_U[self.region] ** 2)) - IQ_reg_err = np.sqrt(np.sum(s_IQ[self.region] ** 2)) - IU_reg_err = np.sqrt(np.sum(s_IU[self.region] ** 2)) - QU_reg_err = np.sqrt(np.sum(s_QU[self.region] ** 2)) + I_reg_err = np.sqrt(np.sum(self.Stokes_cov[0, 0][self.region])) + Q_reg_err = np.sqrt(np.sum(self.Stokes_cov[1, 1][self.region])) + U_reg_err = np.sqrt(np.sum(self.Stokes_cov[2, 2][self.region])) + IQ_reg_err = np.sqrt(np.sum(self.Stokes_cov[0, 1][self.region] ** 2)) + IU_reg_err = np.sqrt(np.sum(self.Stokes_cov[0, 2][self.region] ** 2)) + QU_reg_err = np.sqrt(np.sum(self.Stokes_cov[1, 2][self.region] ** 2)) + I_reg_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[0, 0][self.region])) + Q_reg_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[1, 1][self.region])) + U_reg_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[2, 2][self.region])) + IQ_reg_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[0, 1][self.region] ** 2)) + IU_reg_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[0, 2][self.region] ** 2)) + QU_reg_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[1, 2][self.region] ** 2)) conf = PCconf(QN=Q_reg / I_reg, QN_ERR=Q_reg_err / I_reg, UN=U_reg / I_reg, UN_ERR=U_reg_err / I_reg) if 1.0 - conf > 1e-3: @@ -3432,39 +3463,66 @@ class pol_map(object): ) / I_reg ) + P_reg_stat_err = ( + P_reg + / I_reg + * np.sqrt( + I_reg_stat_err + - 2.0 / (I_reg * P_reg**2) * (Q_reg * IQ_reg_stat_err + U_reg * IU_reg_stat_err) + + 1.0 / (I_reg**2 * P_reg**4) * (Q_reg**2 * Q_reg_stat_err + U_reg**2 * U_reg_stat_err + 2.0 * Q_reg * U_reg * QU_reg_stat_err) + ) + ) + debiased_P_reg = np.sqrt(P_reg**2 - P_reg_stat_err**2) if P_reg**2 > P_reg_stat_err**2 else 0.0 PA_reg = princ_angle((90.0 / np.pi) * np.arctan2(U_reg, Q_reg)) PA_reg_err = (90.0 / (np.pi * (Q_reg**2 + U_reg**2))) * np.sqrt( U_reg**2 * Q_reg_err**2 + Q_reg**2 * U_reg_err**2 - 2.0 * Q_reg * U_reg * QU_reg_err ) - new_cut = np.logical_and(self.region, self.cut) - I_cut = self.I[new_cut].sum() - Q_cut = self.Q[new_cut].sum() - U_cut = self.U[new_cut].sum() - I_cut_err = np.sqrt(np.sum(s_I[new_cut] ** 2)) - Q_cut_err = np.sqrt(np.sum(s_Q[new_cut] ** 2)) - U_cut_err = np.sqrt(np.sum(s_U[new_cut] ** 2)) - IQ_cut_err = np.sqrt(np.sum(s_IQ[new_cut] ** 2)) - IU_cut_err = np.sqrt(np.sum(s_IU[new_cut] ** 2)) - QU_cut_err = np.sqrt(np.sum(s_QU[new_cut] ** 2)) + if cut: + new_cut = np.logical_and(self.region, self.cut) + I_cut = self.I[new_cut].sum() + Q_cut = self.Q[new_cut].sum() + U_cut = self.U[new_cut].sum() + I_cut_err = np.sqrt(np.sum(self.Stokes_cov[0, 0][new_cut])) + Q_cut_err = np.sqrt(np.sum(self.Stokes_cov[1, 1][new_cut])) + U_cut_err = np.sqrt(np.sum(self.Stokes_cov[2, 2][new_cut])) + IQ_cut_err = np.sqrt(np.sum(self.Stokes_cov[0, 1][new_cut] ** 2)) + IU_cut_err = np.sqrt(np.sum(self.Stokes_cov[0, 2][new_cut] ** 2)) + QU_cut_err = np.sqrt(np.sum(self.Stokes_cov[1, 2][new_cut] ** 2)) + I_cut_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[0, 0][new_cut])) + Q_cut_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[1, 1][new_cut])) + U_cut_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[2, 2][new_cut])) + IQ_cut_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[0, 1][new_cut] ** 2)) + IU_cut_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[0, 2][new_cut] ** 2)) + QU_cut_stat_err = np.sqrt(np.sum(self.Stokes_cov_stat[1, 2][new_cut] ** 2)) - with np.errstate(divide="ignore", invalid="ignore"): - P_cut = np.sqrt(Q_cut**2 + U_cut**2) / I_cut - P_cut_err = ( - np.sqrt( - (Q_cut**2 * Q_cut_err**2 + U_cut**2 * U_cut_err**2 + 2.0 * Q_cut * U_cut * QU_cut_err) / (Q_cut**2 + U_cut**2) - + ((Q_cut / I_cut) ** 2 + (U_cut / I_cut) ** 2) * I_cut_err**2 - - 2.0 * (Q_cut / I_cut) * IQ_cut_err - - 2.0 * (U_cut / I_cut) * IU_cut_err + with np.errstate(divide="ignore", invalid="ignore"): + P_cut = np.sqrt(Q_cut**2 + U_cut**2) / I_cut + P_cut_err = ( + np.sqrt( + (Q_cut**2 * Q_cut_err**2 + U_cut**2 * U_cut_err**2 + 2.0 * Q_cut * U_cut * QU_cut_err) / (Q_cut**2 + U_cut**2) + + ((Q_cut / I_cut) ** 2 + (U_cut / I_cut) ** 2) * I_cut_err**2 + - 2.0 * (Q_cut / I_cut) * IQ_cut_err + - 2.0 * (U_cut / I_cut) * IU_cut_err + ) + / I_cut ) - / I_cut - ) + P_cut_stat_err = ( + P_cut + / I_cut + * np.sqrt( + I_cut_stat_err + - 2.0 / (I_cut * P_cut**2) * (Q_cut * IQ_cut_stat_err + U_cut * IU_cut_stat_err) + + 1.0 / (I_cut**2 * P_cut**4) * (Q_cut**2 * Q_cut_stat_err + U_cut**2 * U_cut_stat_err + 2.0 * Q_cut * U_cut * QU_cut_stat_err) + ) + ) + debiased_P_cut = np.sqrt(P_cut**2 - P_cut_stat_err**2) if P_cut**2 > P_cut_stat_err**2 else 0.0 - PA_cut = princ_angle((90.0 / np.pi) * np.arctan2(U_cut, Q_cut)) - PA_cut_err = (90.0 / (np.pi * (Q_cut**2 + U_cut**2))) * np.sqrt( - U_cut**2 * Q_cut_err**2 + Q_cut**2 * U_cut_err**2 - 2.0 * Q_cut * U_cut * QU_cut_err - ) + PA_cut = princ_angle((90.0 / np.pi) * np.arctan2(U_cut, Q_cut)) + PA_cut_err = (90.0 / (np.pi * (Q_cut**2 + U_cut**2))) * np.sqrt( + U_cut**2 * Q_cut_err**2 + Q_cut**2 * U_cut_err**2 - 2.0 * Q_cut * U_cut * QU_cut_err + ) if hasattr(self, "cont"): self.cont.remove() @@ -3480,22 +3538,23 @@ class pol_map(object): self.pivot_wav, sci_not(I_reg * self.map_convert, I_reg_err * self.map_convert, 2) ) + "\n" - + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg * 100.0, np.ceil(P_reg_err * 1000.0) / 10.0) + + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(debiased_P_reg * 100.0, np.ceil(P_reg_err * 1000.0) / 10.0) + "\n" + r"$\Psi^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err * 10.0) / 10.0) + str_conf ) self.str_cut = "" - # self.str_cut = ( - # "\n" - # + r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( - # self.pivot_wav, sci_not(I_cut * self.map_convert, I_cut_err * self.map_convert, 2) - # ) - # + "\n" - # + r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut * 100.0, np.ceil(P_cut_err * 1000.0) / 10.0) - # + "\n" - # + r"$\Psi^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err * 10.0) / 10.0) - # ) + if cut: + self.str_cut = ( + "\n" + + r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( + self.pivot_wav, sci_not(I_cut * self.map_convert, I_cut_err * self.map_convert, 2) + ) + + "\n" + + r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(debiased_P_cut * 100.0, np.ceil(P_cut_err * 1000.0) / 10.0) + + "\n" + + r"$\Psi^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err * 10.0) / 10.0) + ) self.an_int = ax.annotate( self.str_int + self.str_cut, color="white", @@ -3522,16 +3581,17 @@ class pol_map(object): + str_conf ) str_cut = "" - # str_cut = ( - # "\n" - # + r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( - # self.pivot_wav, sci_not(I_cut * self.map_convert, I_cut_err * self.map_convert, 2) - # ) - # + "\n" - # + r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut * 100.0, np.ceil(P_cut_err * 1000.0) / 10.0) - # + "\n" - # + r"$\Psi^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err * 10.0) / 10.0) - # ) + if cut: + str_cut = ( + "\n" + + r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( + self.pivot_wav, sci_not(I_cut * self.map_convert, I_cut_err * self.map_convert, 2) + ) + + "\n" + + r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut * 100.0, np.ceil(P_cut_err * 1000.0) / 10.0) + + "\n" + + r"$\Psi^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err * 10.0) / 10.0) + ) ax.annotate( str_int + str_cut, color="white", diff --git a/package/lib/reduction.py b/package/lib/reduction.py index 749ff5f..59096f9 100755 --- a/package/lib/reduction.py +++ b/package/lib/reduction.py @@ -1301,12 +1301,12 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale # Statistical error: Poisson noise is assumed sigma_flux = np.array([np.sqrt(flux / head["exptime"]) for flux, head in zip(pol_flux, pol_headers)]) - s_IQU_stat = np.zeros(Stokes_cov.shape) + Stokes_cov_stat = np.zeros(Stokes_cov.shape) for i in range(3): - s_IQU_stat[i, i] = np.sum([coeff_stokes[i, k] ** 2 * sigma_flux[k] ** 2 for k in range(len(sigma_flux))], axis=0) + Stokes_cov_stat[i, i] = np.sum([coeff_stokes[i, k] ** 2 * sigma_flux[k] ** 2 for k in range(len(sigma_flux))], axis=0) for j in [k for k in range(3) if k > i]: - s_IQU_stat[i, j] = np.sum([coeff_stokes[i, k] * coeff_stokes[j, k] * sigma_flux[k] ** 2 for k in range(len(sigma_flux))], axis=0) - s_IQU_stat[j, i] = np.sum([coeff_stokes[i, k] * coeff_stokes[j, k] * sigma_flux[k] ** 2 for k in range(len(sigma_flux))], axis=0) + Stokes_cov_stat[i, j] = np.sum([coeff_stokes[i, k] * coeff_stokes[j, k] * sigma_flux[k] ** 2 for k in range(len(sigma_flux))], axis=0) + Stokes_cov_stat[j, i] = np.sum([coeff_stokes[i, k] * coeff_stokes[j, k] * sigma_flux[k] ** 2 for k in range(len(sigma_flux))], axis=0) # Compute the derivative of each Stokes parameter with respect to the polarizer orientation dIQU_dtheta = np.zeros(Stokes_cov.shape) @@ -1359,19 +1359,19 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale ) # Compute the uncertainty associated with the polarizers' orientation (see Kishimoto 1999) - s_IQU_axis = np.zeros(Stokes_cov.shape) + Stokes_cov_axis = np.zeros(Stokes_cov.shape) for i in range(3): - s_IQU_axis[i, i] = np.sum([dIQU_dtheta[i, k] ** 2 * globals()["sigma_theta"][k] ** 2 for k in range(len(globals()["sigma_theta"]))], axis=0) + Stokes_cov_axis[i, i] = np.sum([dIQU_dtheta[i, k] ** 2 * globals()["sigma_theta"][k] ** 2 for k in range(len(globals()["sigma_theta"]))], axis=0) for j in [k for k in range(3) if k > i]: - s_IQU_axis[i, j] = np.sum( + Stokes_cov_axis[i, j] = np.sum( [dIQU_dtheta[i, k] * dIQU_dtheta[j, k] * globals()["sigma_theta"][k] ** 2 for k in range(len(globals()["sigma_theta"]))], axis=0 ) - s_IQU_axis[j, i] = np.sum( + Stokes_cov_axis[j, i] = np.sum( [dIQU_dtheta[i, k] * dIQU_dtheta[j, k] * globals()["sigma_theta"][k] ** 2 for k in range(len(globals()["sigma_theta"]))], axis=0 ) # Add quadratically the uncertainty to the Stokes covariance matrix - Stokes_cov += s_IQU_axis + s_IQU_stat + Stokes_cov += Stokes_cov_axis + Stokes_cov_stat # Save values to single header header_stokes = pol_headers[0] @@ -1445,10 +1445,10 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale header_stokes["PA_int"] = (PA_diluted, "Integrated polarization angle") header_stokes["sPA_int"] = (np.ceil(PA_diluted_err * 10.0) / 10.0, "Integrated polarization angle error") - return Stokes, Stokes_cov, header_stokes, s_IQU_stat + return Stokes, Stokes_cov, header_stokes, Stokes_cov_stat -def compute_pol(Stokes, Stokes_cov, header_stokes, s_IQU_stat=None): +def compute_pol(Stokes, Stokes_cov, header_stokes, Stokes_cov_stat=None): """ Compute the polarization degree (in %) and angle (in deg) and their respective errors from given Stokes parameters. @@ -1524,7 +1524,7 @@ def compute_pol(Stokes, Stokes_cov, header_stokes, s_IQU_stat=None): s_P_P = np.ones(Stokes[0].shape) * fmax s_PA_P = np.ones(Stokes[0].shape) * fmax maskP = np.logical_and(mask, P > 0.0) - if s_IQU_stat is not None: + if Stokes_cov_stat is not None: # If IQU covariance matrix containing only statistical error is given propagate to P and PA # Catch Invalid value in sqrt when diagonal terms are big with warnings.catch_warnings(record=True) as _: @@ -1532,13 +1532,15 @@ def compute_pol(Stokes, Stokes_cov, header_stokes, s_IQU_stat=None): P[maskP] / Stokes[0][maskP] * np.sqrt( - s_IQU_stat[0, 0][maskP] - - 2.0 / (Stokes[0][maskP] * P[maskP] ** 2) * (Stokes[1][maskP] * s_IQU_stat[0, 1][maskP] + Stokes[2][maskP] * s_IQU_stat[0, 2][maskP]) + Stokes_cov_stat[0, 0][maskP] + - 2.0 + / (Stokes[0][maskP] * P[maskP] ** 2) + * (Stokes[1][maskP] * Stokes_cov_stat[0, 1][maskP] + Stokes[2][maskP] * Stokes_cov_stat[0, 2][maskP]) + 1.0 / (Stokes[0][maskP] ** 2 * P[maskP] ** 4) * ( - Stokes[1][maskP] ** 2 * s_IQU_stat[1, 1][maskP] - + Stokes[2][maskP] ** 2 * s_IQU_stat[2, 2][maskP] * Stokes[1][maskP] * Stokes[2][maskP] * s_IQU_stat[1, 2][maskP] + Stokes[1][maskP] ** 2 * Stokes_cov_stat[1, 1][maskP] + + Stokes[2][maskP] ** 2 * Stokes_cov_stat[2, 2][maskP] * Stokes[1][maskP] * Stokes[2][maskP] * Stokes_cov_stat[1, 2][maskP] ) ) ) @@ -1546,9 +1548,9 @@ def compute_pol(Stokes, Stokes_cov, header_stokes, s_IQU_stat=None): 90.0 / (np.pi * Stokes[0][maskP] ** 2 * P[maskP] ** 2) * ( - Stokes[1][maskP] ** 2 * s_IQU_stat[2, 2][maskP] - + Stokes[2][maskP] * s_IQU_stat[1, 1][maskP] - - 2.0 * Stokes[1][maskP] * Stokes[2][maskP] * s_IQU_stat[1, 2][maskP] + Stokes[1][maskP] ** 2 * Stokes_cov_stat[2, 2][maskP] + + Stokes[2][maskP] * Stokes_cov_stat[1, 1][maskP] + - 2.0 * Stokes[1][maskP] * Stokes[2][maskP] * Stokes_cov_stat[1, 2][maskP] ) ) else: @@ -1576,7 +1578,7 @@ def compute_pol(Stokes, Stokes_cov, header_stokes, s_IQU_stat=None): return P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P -def rotate_Stokes(Stokes, Stokes_cov, data_mask, header_stokes, s_IQU_stat=None, SNRi_cut=None): +def rotate_Stokes(Stokes, Stokes_cov, data_mask, header_stokes, Stokes_cov_stat=None, SNRi_cut=None): """ Use scipy.ndimage.rotate to rotate I_stokes to an angle, and a rotation matrix to rotate Q, U of a given angle in degrees and update header @@ -1642,16 +1644,16 @@ def rotate_Stokes(Stokes, Stokes_cov, data_mask, header_stokes, s_IQU_stat=None, new_Stokes[:3, i, j] = np.dot(mrot, new_Stokes[:3, i, j]) new_Stokes_cov[:3, :3, i, j] = np.dot(mrot, np.dot(new_Stokes_cov[:3, :3, i, j], mrot.T)) - if s_IQU_stat is not None: - s_IQU_stat = zeropad(s_IQU_stat, [*s_IQU_stat.shape[:-2], *shape]) - new_s_IQU_stat = np.zeros((*s_IQU_stat.shape[:-2], *shape)) + if Stokes_cov_stat is not None: + Stokes_cov_stat = zeropad(Stokes_cov_stat, [*Stokes_cov_stat.shape[:-2], *shape]) + new_Stokes_cov_stat = np.zeros((*Stokes_cov_stat.shape[:-2], *shape)) for i in range(3): for j in range(3): - new_s_IQU_stat[i, j] = sc_rotate(s_IQU_stat[i, j], ang, order=1, reshape=False, cval=0.0) - new_s_IQU_stat[i, i] = np.abs(new_s_IQU_stat[i, i]) + new_Stokes_cov_stat[i, j] = sc_rotate(Stokes_cov_stat[i, j], ang, order=1, reshape=False, cval=0.0) + new_Stokes_cov_stat[i, i] = np.abs(new_Stokes_cov_stat[i, i]) for i in range(shape[0]): for j in range(shape[1]): - new_s_IQU_stat[:3, :3, i, j] = np.dot(mrot, np.dot(new_s_IQU_stat[:3, :3, i, j], mrot.T)) + new_Stokes_cov_stat[:3, :3, i, j] = np.dot(mrot, np.dot(new_Stokes_cov_stat[:3, :3, i, j], mrot.T)) # Update headers to new angle mrot = np.array([[np.cos(-alpha), -np.sin(-alpha)], [np.sin(-alpha), np.cos(-alpha)]]) @@ -1701,8 +1703,8 @@ def rotate_Stokes(Stokes, Stokes_cov, data_mask, header_stokes, s_IQU_stat=None, new_header_stokes["PA_int"] = (PA_diluted, "Integrated polarization angle") new_header_stokes["sPA_int"] = (np.ceil(PA_diluted_err * 10.0) / 10.0, "Integrated polarization angle error") - if s_IQU_stat is not None: - return new_Stokes, new_Stokes_cov, new_data_mask, new_header_stokes, new_s_IQU_stat + if Stokes_cov_stat is not None: + return new_Stokes, new_Stokes_cov, new_data_mask, new_header_stokes, new_Stokes_cov_stat else: return new_Stokes, new_Stokes_cov, new_data_mask, new_header_stokes