diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index c4d7272..37b52c2 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -40,13 +40,13 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= display_crop = False # Background estimation - error_sub_type = "scott" # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) - subtract_error = 0.50 + error_sub_type = "freedman-diaconis" # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) + subtract_error = 1.33 display_bkg = True # Data binning - pxsize = 4 - pxscale = "px" # pixel, arcsec or full + pxsize = 0.05 + pxscale = "arcsec" # pixel, arcsec or full rebin_operation = "sum" # sum or average # Alignement @@ -59,8 +59,8 @@ 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 = 1.5 # If None, no smoothing is done - smoothing_scale = "px" # pixel or arcsec + smoothing_FWHM = 0.075 # If None, no smoothing is done + smoothing_scale = "arcsec" # pixel or arcsec # Rotation rotate_North = True @@ -216,28 +216,26 @@ 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 - I_stokes, Q_stokes, U_stokes, Stokes_cov, header_stokes, s_IQU_stat = proj_red.compute_Stokes( + I_stokes, Q_stokes, U_stokes, Stokes_cov, Stokes_stat_cov, header_stokes = proj_red.compute_Stokes( 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, s_IQU_stat_bkg = proj_red.compute_Stokes( + I_bkg, Q_bkg, U_bkg, S_cov_bkg, S_stat_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 3: # Rotate images to have North up if rotate_North: - I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_stokes, s_IQU_stat = proj_red.rotate_Stokes( - I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_stokes, s_IQU_stat=s_IQU_stat, SNRi_cut=None + I_stokes, Q_stokes, U_stokes, Stokes_cov, Stokes_stat_cov, data_mask, header_stokes = proj_red.rotate_Stokes( + I_stokes, Q_stokes, U_stokes, Stokes_cov, Stokes_stat_cov, data_mask, header_stokes, SNRi_cut=None ) - I_bkg, Q_bkg, U_bkg, S_cov_bkg, data_mask_bkg, header_bkg, s_IQU_stat_bkg = proj_red.rotate_Stokes( - I_bkg, Q_bkg, U_bkg, S_cov_bkg, np.array(True).reshape(1, 1), header_bkg, s_IQU_stat=s_IQU_stat_bkg, SNRi_cut=None + I_bkg, Q_bkg, U_bkg, S_cov_bkg, S_stat_cov_bkg, data_mask_bkg, header_bkg = proj_red.rotate_Stokes( + I_bkg, Q_bkg, U_bkg, S_cov_bkg, S_stat_cov_bkg, np.array(True).reshape(1, 1), header_bkg, 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(I_stokes, Q_stokes, U_stokes, Stokes_cov, header_stokes, s_IQU_stat=s_IQU_stat) - 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, s_IQU_stat=s_IQU_stat_bkg - ) + 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, Stokes_stat_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, S_stat_cov_bkg, header_bkg) # Step 4: # Save image to FITS. @@ -247,6 +245,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= Q_stokes, U_stokes, Stokes_cov, + Stokes_stat_cov, P, debiased_P, s_P, diff --git a/package/lib/fits.py b/package/lib/fits.py index 43d62d9..898fa77 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -106,7 +106,23 @@ 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, + Stokes_stat_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, ): """ Save computed polarimetry parameters to a single fits file, @@ -186,11 +202,15 @@ def save_Stokes( s_PA_P = s_PA_P[vertex[2] : vertex[3], vertex[0] : vertex[1]] new_Stokes_cov = np.zeros((*Stokes_cov.shape[:-2], *shape[::-1])) + new_Stokes_stat_cov = np.zeros((*Stokes_stat_cov.shape[:-2], *shape[::-1])) for i in range(3): for j in range(3): Stokes_cov[i, j][(1 - data_mask).astype(bool)] = 0.0 new_Stokes_cov[i, j] = Stokes_cov[i, j][vertex[2] : vertex[3], vertex[0] : vertex[1]] + Stokes_stat_cov[i, j][(1 - data_mask).astype(bool)] = 0.0 + new_Stokes_stat_cov[i, j] = Stokes_stat_cov[i, j][vertex[2] : vertex[3], vertex[0] : vertex[1]] Stokes_cov = new_Stokes_cov + Stokes_stat_cov = new_Stokes_stat_cov data_mask = data_mask[vertex[2] : vertex[3], vertex[0] : vertex[1]] data_mask = data_mask.astype(float, copy=False) @@ -210,6 +230,7 @@ def save_Stokes( [Q_stokes, "Q_stokes"], [U_stokes, "U_stokes"], [Stokes_cov, "IQU_cov_matrix"], + [Stokes_stat_cov, "IQU_stat_cov_matrix"], [P, "Pol_deg"], [debiased_P, "Pol_deg_debiased"], [s_P, "Pol_deg_err"], @@ -221,7 +242,7 @@ def save_Stokes( ]: hdu_header = header.copy() hdu_header["datatype"] = name - if not name == "IQU_cov_matrix": + if not name[-10:] == "cov_matrix": data[(1 - data_mask).astype(bool)] = 0.0 hdu = fits.ImageHDU(data=data, header=hdu_header) hdu.name = name diff --git a/package/lib/plots.py b/package/lib/plots.py index 5d1565a..c614263 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -2079,7 +2079,7 @@ class crop_Stokes(crop_map): # Crop dataset for dataset in self.hdul_crop: - if dataset.header["datatype"] == "IQU_cov_matrix": + if dataset.header["datatype"][-10:] == "cov_matrix": stokes_cov = np.zeros((3, 3, shape[1], shape[0])) for i in range(3): for j in range(3): @@ -2104,16 +2104,22 @@ class crop_Stokes(crop_map): if self.fig.canvas.manager.toolbar.mode == "": self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, button=[1], spancoords="pixels", useblit=True) # Update integrated values - mask = np.logical_and(self.hdul_crop["data_mask"].data.astype(bool), self.hdul_crop[0].data > 0) - I_diluted = self.hdul_crop["i_stokes"].data[mask].sum() - Q_diluted = self.hdul_crop["q_stokes"].data[mask].sum() - U_diluted = self.hdul_crop["u_stokes"].data[mask].sum() - I_diluted_err = np.sqrt(np.sum(self.hdul_crop["iqu_cov_matrix"].data[0, 0][mask])) - Q_diluted_err = np.sqrt(np.sum(self.hdul_crop["iqu_cov_matrix"].data[1, 1][mask])) - U_diluted_err = np.sqrt(np.sum(self.hdul_crop["iqu_cov_matrix"].data[2, 2][mask])) - IQ_diluted_err = np.sqrt(np.sum(self.hdul_crop["iqu_cov_matrix"].data[0, 1][mask] ** 2)) - IU_diluted_err = np.sqrt(np.sum(self.hdul_crop["iqu_cov_matrix"].data[0, 2][mask] ** 2)) - QU_diluted_err = np.sqrt(np.sum(self.hdul_crop["iqu_cov_matrix"].data[1, 2][mask] ** 2)) + mask = np.logical_and(self.hdul_crop["DATA_MASK"].data.astype(bool), self.hdul_crop[0].data > 0) + I_diluted = self.hdul_crop["I_STOKES"].data[mask].sum() + Q_diluted = self.hdul_crop["Q_STOKES"].data[mask].sum() + U_diluted = self.hdul_crop["U_STOKES"].data[mask].sum() + I_diluted_err = np.sqrt(np.sum(self.hdul_crop["IQU_COV_MATRIX"].data[0, 0][mask])) + Q_diluted_err = np.sqrt(np.sum(self.hdul_crop["IQU_COV_MATRIX"].data[1, 1][mask])) + U_diluted_err = np.sqrt(np.sum(self.hdul_crop["IQU_COV_MATRIX"].data[2, 2][mask])) + IQ_diluted_err = np.sqrt(np.sum(self.hdul_crop["IQU_COV_MATRIX"].data[0, 1][mask] ** 2)) + IU_diluted_err = np.sqrt(np.sum(self.hdul_crop["IQU_COV_MATRIX"].data[0, 2][mask] ** 2)) + QU_diluted_err = np.sqrt(np.sum(self.hdul_crop["IQU_COV_MATRIX"].data[1, 2][mask] ** 2)) + I_diluted_stat_err = np.sqrt(np.sum(self.hdul_crop["IQU_STAT_COV_MATRIX"].data[0, 0][mask])) + Q_diluted_stat_err = np.sqrt(np.sum(self.hdul_crop["IQU_STAT_COV_MATRIX"].data[1, 1][mask])) + U_diluted_stat_err = np.sqrt(np.sum(self.hdul_crop["IQU_STAT_COV_MATRIX"].data[2, 2][mask])) + IQ_diluted_stat_err = np.sqrt(np.sum(self.hdul_crop["IQU_STAT_COV_MATRIX"].data[0, 1][mask] ** 2)) + IU_diluted_stat_err = np.sqrt(np.sum(self.hdul_crop["IQU_STAT_COV_MATRIX"].data[0, 2][mask] ** 2)) + QU_diluted_stat_err = np.sqrt(np.sum(self.hdul_crop["IQU_STAT_COV_MATRIX"].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( @@ -2122,6 +2128,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( @@ -2131,7 +2149,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") @@ -2961,6 +2979,10 @@ class pol_map(object): def IQU_cov(self): return self.Stokes["IQU_COV_MATRIX"].data + @property + def IQU_stat_cov(self): + return self.Stokes["IQU_STAT_COV_MATRIX"].data + @property def P(self): return self.Stokes["POL_DEG_DEBIASED"].data @@ -3283,27 +3305,26 @@ class pol_map(object): s_I = np.sqrt(self.IQU_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] - 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)) + I_cut_err = np.sqrt(np.sum(self.IQU_cov[0, 0][self.cut])) + Q_cut_err = np.sqrt(np.sum(self.IQU_cov[1, 1][self.cut])) + U_cut_err = np.sqrt(np.sum(self.IQU_cov[2, 2][self.cut])) + IQ_cut_err = np.sqrt(np.sum(self.IQU_cov[0, 1][self.cut] ** 2)) + IU_cut_err = np.sqrt(np.sum(self.IQU_cov[0, 2][self.cut] ** 2)) + QU_cut_err = np.sqrt(np.sum(self.IQU_cov[1, 2][self.cut] ** 2)) + I_cut_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[0, 0][self.cut])) + Q_cut_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[1, 1][self.cut])) + U_cut_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[2, 2][self.cut])) + IQ_cut_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[0, 1][self.cut] ** 2)) + IU_cut_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[0, 2][self.cut] ** 2)) + QU_cut_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[1, 2][self.cut] ** 2)) with np.errstate(divide="ignore", invalid="ignore"): P_cut = np.sqrt(Q_cut**2 + U_cut**2) / I_cut @@ -3316,6 +3337,18 @@ class pol_map(object): ) / 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) + ) + ) + mask = P_cut**2 > P_cut_stat_err + debiased_P_cut = np.zeros(P_cut.shape) + debiased_P_cut[mask] = np.sqrt(P_cut[mask] ** 2 - P_cut_stat_err[mask] ** 2) 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( @@ -3323,22 +3356,21 @@ class pol_map(object): ) 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.IQU_cov[0, 0][self.region])) + Q_reg_err = np.sqrt(np.sum(self.IQU_cov[1, 1][self.region])) + U_reg_err = np.sqrt(np.sum(self.IQU_cov[2, 2][self.region])) + IQ_reg_err = np.sqrt(np.sum(self.IQU_cov[0, 1][self.region] ** 2)) + IU_reg_err = np.sqrt(np.sum(self.IQU_cov[0, 2][self.region] ** 2)) + QU_reg_err = np.sqrt(np.sum(self.IQU_cov[1, 2][self.region] ** 2)) + I_reg_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[0, 0][self.region])) + Q_reg_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[1, 1][self.region])) + U_reg_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[2, 2][self.region])) + IQ_reg_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[0, 1][self.region] ** 2)) + IU_reg_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[0, 2][self.region] ** 2)) + QU_reg_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[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: @@ -3355,6 +3387,16 @@ 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( @@ -3365,12 +3407,18 @@ class pol_map(object): 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)) + I_cut_err = np.sqrt(np.sum(self.IQU_cov[0, 0][new_cut])) + Q_cut_err = np.sqrt(np.sum(self.IQU_cov[1, 1][new_cut])) + U_cut_err = np.sqrt(np.sum(self.IQU_cov[2, 2][new_cut])) + IQ_cut_err = np.sqrt(np.sum(self.IQU_cov[0, 1][new_cut] ** 2)) + IU_cut_err = np.sqrt(np.sum(self.IQU_cov[0, 2][new_cut] ** 2)) + QU_cut_err = np.sqrt(np.sum(self.IQU_cov[1, 2][new_cut] ** 2)) + I_cut_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[0, 0][new_cut])) + Q_cut_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[1, 1][new_cut])) + U_cut_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[2, 2][new_cut])) + IQ_cut_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[0, 1][new_cut] ** 2)) + IU_cut_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[0, 2][new_cut] ** 2)) + QU_cut_stat_err = np.sqrt(np.sum(self.IQU_stat_cov[1, 2][new_cut] ** 2)) with np.errstate(divide="ignore", invalid="ignore"): P_cut = np.sqrt(Q_cut**2 + U_cut**2) / I_cut @@ -3383,6 +3431,18 @@ class pol_map(object): ) / 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) + ) + ) + mask = P_cut**2 > P_cut_stat_err + debiased_P_cut = np.zeros(P_cut.shape) + debiased_P_cut[mask] = np.sqrt(P_cut[mask] ** 2 - P_cut_stat_err[mask] ** 2) 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( @@ -3403,7 +3463,7 @@ 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"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err * 10.0) / 10.0) + str_conf @@ -3415,7 +3475,7 @@ class pol_map(object): # 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) + # + 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"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err * 10.0) / 10.0) # ) @@ -3439,7 +3499,7 @@ 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"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg, np.ceil(PA_reg_err * 10.0) / 10.0) + str_conf @@ -3451,7 +3511,7 @@ class pol_map(object): # 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) + # + 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"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut, np.ceil(PA_cut_err * 10.0) / 10.0) # ) diff --git a/package/lib/reduction.py b/package/lib/reduction.py index 24f788f..6fe8ecd 100755 --- a/package/lib/reduction.py +++ b/package/lib/reduction.py @@ -1310,12 +1310,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_stat_cov = np.zeros(Stokes_cov.shape) for i in range(Stokes_cov.shape[0]): - 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_stat_cov[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([abs(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([abs(coeff_stokes[i, k] * coeff_stokes[j, k]) * sigma_flux[k] ** 2 for k in range(len(sigma_flux))], axis=0) + Stokes_stat_cov[i, j] = np.sum([abs(coeff_stokes[i, k] * coeff_stokes[j, k]) * sigma_flux[k] ** 2 for k in range(len(sigma_flux))], axis=0) + Stokes_stat_cov[j, i] = np.sum([abs(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) @@ -1368,19 +1368,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_axis_cov = np.zeros(Stokes_cov.shape) for i in range(Stokes_cov.shape[0]): - 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_axis_cov[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_axis_cov[i, j] = np.sum( [abs(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_axis_cov[j, i] = np.sum( [abs(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_axis_cov + Stokes_stat_cov # Save values to single header header_stokes = pol_headers[0] @@ -1414,8 +1414,8 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale for i in range(3): Stokes_cov[i, i] = np.sum([exp**2 * cov for exp, cov in zip(all_exp, all_Stokes_cov[:, i, i])], axis=0) / all_exp.sum() ** 2 for j in [x for x in range(3) if x != i]: - Stokes_cov[i, j] = np.sqrt(np.sum([exp**2 * cov**2 for exp, cov in zip(all_exp, all_Stokes_cov[:, i, j])], axis=0) / all_exp.sum() ** 2) - Stokes_cov[j, i] = np.sqrt(np.sum([exp**2 * cov**2 for exp, cov in zip(all_exp, all_Stokes_cov[:, j, i])], axis=0) / all_exp.sum() ** 2) + Stokes_cov[i, j] = np.sum([exp**2 * cov**2 for exp, cov in zip(all_exp, all_Stokes_cov[:, i, j])], axis=0) / all_exp.sum() ** 2 + Stokes_cov[j, i] = np.sum([exp**2 * cov**2 for exp, cov in zip(all_exp, all_Stokes_cov[:, j, i])], axis=0) / all_exp.sum() ** 2 # Save values to single header header_stokes = all_header_stokes[0] @@ -1430,6 +1430,7 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale Q_stokes[np.isnan(Q_stokes)] = 0.0 U_stokes[np.isnan(U_stokes)] = 0.0 Stokes_cov[np.isnan(Stokes_cov)] = fmax + Stokes_stat_cov[np.isnan(Stokes_cov)] = fmax if integrate: # Compute integrated values for P, PA before any rotation @@ -1443,6 +1444,12 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale IQ_diluted_err = np.sqrt(np.sum(Stokes_cov[0, 1][mask] ** 2)) IU_diluted_err = np.sqrt(np.sum(Stokes_cov[0, 2][mask] ** 2)) QU_diluted_err = np.sqrt(np.sum(Stokes_cov[1, 2][mask] ** 2)) + I_diluted_stat_err = np.sqrt(np.sum(Stokes_stat_cov[0, 0][mask])) + Q_diluted_stat_err = np.sqrt(np.sum(Stokes_stat_cov[1, 1][mask])) + U_diluted_stat_err = np.sqrt(np.sum(Stokes_stat_cov[2, 2][mask])) + IQ_diluted_stat_err = np.sqrt(np.sum(Stokes_stat_cov[0, 1][mask] ** 2)) + IU_diluted_stat_err = np.sqrt(np.sum(Stokes_stat_cov[0, 2][mask] ** 2)) + QU_diluted_stat_err = np.sqrt(np.sum(Stokes_stat_cov[1, 2][mask] ** 2)) P_diluted = np.sqrt(Q_diluted**2 + U_diluted**2) / I_diluted P_diluted_err = np.sqrt( @@ -1451,21 +1458,33 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale - 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( U_diluted**2 * Q_diluted_err**2 + Q_diluted**2 * U_diluted_err**2 - 2.0 * Q_diluted * U_diluted * QU_diluted_err ) - header_stokes["P_int"] = (P_diluted, "Integrated polarization degree") + header_stokes["P_int"] = (debiased_P_diluted, "Integrated polarization degree") header_stokes["sP_int"] = (np.ceil(P_diluted_err * 1000.0) / 1000.0, "Integrated polarization degree error") 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 I_stokes, Q_stokes, U_stokes, Stokes_cov, header_stokes, s_IQU_stat + return I_stokes, Q_stokes, U_stokes, Stokes_cov, Stokes_stat_cov, header_stokes -def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, header_stokes, s_IQU_stat=None): +def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, Stokes_stat_cov, header_stokes): """ Compute the polarization degree (in %) and angle (in deg) and their respective errors from given Stokes parameters. @@ -1540,45 +1559,34 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, header_stokes, s_IQU_s s_P[np.isnan(s_P)] = fmax s_PA[np.isnan(s_PA)] = fmax - # Compute the total exposure time so that - # I_stokes*exp_tot = N_tot the total number of events - N_obs = I_stokes * float(header_stokes["exptime"]) - # Errors on P, PA supposing Poisson noise s_P_P = np.ones(I_stokes.shape) * fmax s_PA_P = np.ones(I_stokes.shape) * fmax maskP = np.logical_and(mask, P > 0.0) - if s_IQU_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 - s_P_P[maskP] = ( - P[maskP] - / I_stokes[maskP] - * np.sqrt( - s_IQU_stat[0, 0][maskP] - - 2.0 / (I_stokes[maskP] * P[maskP] ** 2) * (Q_stokes[maskP] * s_IQU_stat[0, 1][maskP] + U_stokes[maskP] * s_IQU_stat[0, 2][maskP]) - + 1.0 - / (I_stokes[maskP] ** 2 * P[maskP] ** 4) - * ( - Q_stokes[maskP] ** 2 * s_IQU_stat[1, 1][maskP] - + U_stokes[maskP] ** 2 * s_IQU_stat[2, 2][maskP] - + 2.0 * Q_stokes[maskP] * U_stokes[maskP] * s_IQU_stat[1, 2][maskP] - ) - ) - ) - s_PA_P[maskP] = ( - 90.0 - / (np.pi * I_stokes[maskP] ** 2 * P[maskP] ** 2) + s_P_P[maskP] = ( + P[maskP] + / I_stokes[maskP] + * np.sqrt( + Stokes_stat_cov[0, 0][maskP] + - 2.0 / (I_stokes[maskP] * P[maskP] ** 2) * (Q_stokes[maskP] * Stokes_stat_cov[0, 1][maskP] + U_stokes[maskP] * Stokes_stat_cov[0, 2][maskP]) + + 1.0 + / (I_stokes[maskP] ** 2 * P[maskP] ** 4) * ( - Q_stokes[maskP] ** 2 * s_IQU_stat[2, 2][maskP] - + U_stokes[maskP] * s_IQU_stat[1, 1][maskP] - - 2.0 * Q_stokes[maskP] * U_stokes[maskP] * s_IQU_stat[1, 2][maskP] + Q_stokes[maskP] ** 2 * Stokes_stat_cov[1, 1][maskP] + + U_stokes[maskP] ** 2 * Stokes_stat_cov[2, 2][maskP] + + 2.0 * Q_stokes[maskP] * U_stokes[maskP] * Stokes_stat_cov[1, 2][maskP] ) ) - else: - # Approximate Poisson error for P and PA - s_P_P[mask] = np.sqrt(2.0 / N_obs[mask]) - s_PA_P[maskP] = s_P_P[maskP] / P[maskP] * 90.0 / np.pi + ) + s_PA_P[maskP] = ( + 90.0 + / (np.pi * I_stokes[maskP] ** 2 * P[maskP] ** 2) + * ( + Q_stokes[maskP] ** 2 * Stokes_stat_cov[2, 2][maskP] + + U_stokes[maskP] * Stokes_stat_cov[1, 1][maskP] + - 2.0 * Q_stokes[maskP] * U_stokes[maskP] * Stokes_stat_cov[1, 2][maskP] + ) + ) # Catch expected "OverflowWarning" as wrong pixel have an overflowing error with warnings.catch_warnings(record=True) as _: @@ -1600,7 +1608,7 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, header_stokes, s_IQU_s return P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P -def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_stokes, s_IQU_stat=None, SNRi_cut=None): +def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, Stokes_stat_cov, data_mask, header_stokes, 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 @@ -1617,7 +1625,11 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_st Image (2D floats) containing the Stokes parameters accounting for +45/-45deg linear polarization intensity Stokes_cov : numpy.ndarray - Covariance matrix of the Stokes parameters I, Q, U. + Covariance matrix containing all uncertainties of the Stokes + parameters I, Q, U. + Stokes_stat_cov : numpy.ndarray + Covariance matrix containing statistical uncertainty of the Stokes + parameters I, Q, U. data_mask : numpy.ndarray 2D boolean array delimiting the data to work on. header_stokes : astropy.fits.header.Header @@ -1639,6 +1651,8 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_st accounting for +45/-45deg linear polarization intensity. new_Stokes_cov : numpy.ndarray Updated covariance matrix of the Stokes parameters I, Q, U. + new_Stokes_stat_cov : numpy.ndarray + Updated statistical covariance matrix of the Stokes parameters I, Q, U. new_header_stokes : astropy.fits.header.Header Updated Header file associated with the Stokes fluxes accounting for the new orientation angle. @@ -1670,11 +1684,9 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_st Q_stokes = zeropad(Q_stokes, shape) U_stokes = zeropad(U_stokes, shape) data_mask = zeropad(data_mask, shape) - Stokes_cov = zeropad(Stokes_cov, [*Stokes_cov.shape[:-2], *shape]) new_I_stokes = np.zeros(shape) new_Q_stokes = np.zeros(shape) new_U_stokes = np.zeros(shape) - new_Stokes_cov = np.zeros((*Stokes_cov.shape[:-2], *shape)) # Rotate original images using scipy.ndimage.rotate new_I_stokes = sc_rotate(I_stokes, ang, order=1, reshape=False, cval=0.0) @@ -1683,6 +1695,10 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_st new_data_mask = sc_rotate(data_mask.astype(float) * 10.0, ang, order=1, reshape=False, cval=0.0) new_data_mask[new_data_mask < 1.0] = 0.0 new_data_mask = new_data_mask.astype(bool) + + # Rotate covariance matrix + Stokes_cov = zeropad(Stokes_cov, [*Stokes_cov.shape[:-2], *shape]) + new_Stokes_cov = np.zeros((*Stokes_cov.shape[:-2], *shape)) for i in range(3): for j in range(3): new_Stokes_cov[i, j] = sc_rotate(Stokes_cov[i, j], ang, order=1, reshape=False, cval=0.0) @@ -1693,16 +1709,16 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_st new_I_stokes[i, j], new_Q_stokes[i, j], new_U_stokes[i, j] = np.dot(mrot, np.array([new_I_stokes[i, j], new_Q_stokes[i, j], new_U_stokes[i, j]])).T new_Stokes_cov[:, :, i, j] = np.dot(mrot, np.dot(new_Stokes_cov[:, :, 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)) - 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]) - for i in range(shape[0]): - for j in range(shape[1]): - new_s_IQU_stat[:, :, i, j] = np.dot(mrot, np.dot(new_s_IQU_stat[:, :, i, j], mrot.T)) + # Rotate statistical covariance matrix + Stokes_stat_cov = zeropad(Stokes_stat_cov, [*Stokes_stat_cov.shape[:-2], *shape]) + new_Stokes_stat_cov = np.zeros((*Stokes_stat_cov.shape[:-2], *shape)) + for i in range(3): + for j in range(3): + new_Stokes_stat_cov[i, j] = sc_rotate(Stokes_stat_cov[i, j], ang, order=1, reshape=False, cval=0.0) + new_Stokes_stat_cov[i, i] = np.abs(new_Stokes_stat_cov[i, i]) + for i in range(shape[0]): + for j in range(shape[1]): + new_Stokes_stat_cov[:, :, i, j] = np.dot(mrot, np.dot(new_Stokes_stat_cov[:, :, i, j], mrot.T)) # Update headers to new angle mrot = np.array([[np.cos(-alpha), -np.sin(-alpha)], [np.sin(-alpha), np.cos(-alpha)]]) @@ -1732,35 +1748,50 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, header_st I_diluted = new_I_stokes[mask].sum() Q_diluted = new_Q_stokes[mask].sum() U_diluted = new_U_stokes[mask].sum() - I_diluted_err = np.sqrt(np.sum(new_Stokes_cov[0, 0][mask])) - Q_diluted_err = np.sqrt(np.sum(new_Stokes_cov[1, 1][mask])) - U_diluted_err = np.sqrt(np.sum(new_Stokes_cov[2, 2][mask])) - IQ_diluted_err = np.sqrt(np.sum(new_Stokes_cov[0, 1][mask] ** 2)) - IU_diluted_err = np.sqrt(np.sum(new_Stokes_cov[0, 2][mask] ** 2)) - QU_diluted_err = np.sqrt(np.sum(new_Stokes_cov[1, 2][mask] ** 2)) + I_diluted_err = np.sqrt(np.sum(Stokes_cov[0, 0][mask])) + Q_diluted_err = np.sqrt(np.sum(Stokes_cov[1, 1][mask])) + U_diluted_err = np.sqrt(np.sum(Stokes_cov[2, 2][mask])) + IQ_diluted_err = np.sqrt(np.sum(Stokes_cov[0, 1][mask] ** 2)) + IU_diluted_err = np.sqrt(np.sum(Stokes_cov[0, 2][mask] ** 2)) + QU_diluted_err = np.sqrt(np.sum(Stokes_cov[1, 2][mask] ** 2)) + I_diluted_stat_err = np.sqrt(np.sum(Stokes_stat_cov[0, 0][mask])) + Q_diluted_stat_err = np.sqrt(np.sum(Stokes_stat_cov[1, 1][mask])) + U_diluted_stat_err = np.sqrt(np.sum(Stokes_stat_cov[2, 2][mask])) + IQ_diluted_stat_err = np.sqrt(np.sum(Stokes_stat_cov[0, 1][mask] ** 2)) + IU_diluted_stat_err = np.sqrt(np.sum(Stokes_stat_cov[0, 2][mask] ** 2)) + QU_diluted_stat_err = np.sqrt(np.sum(Stokes_stat_cov[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( + P_diluted_err = np.sqrt( (Q_diluted**2 * Q_diluted_err**2 + U_diluted**2 * U_diluted_err**2 + 2.0 * Q_diluted * U_diluted * QU_diluted_err) / (Q_diluted**2 + U_diluted**2) + ((Q_diluted / I_diluted) ** 2 + (U_diluted / I_diluted) ** 2) * I_diluted_err**2 - 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( U_diluted**2 * Q_diluted_err**2 + Q_diluted**2 * U_diluted_err**2 - 2.0 * Q_diluted * U_diluted * QU_diluted_err ) - new_header_stokes["P_int"] = (P_diluted, "Integrated polarization degree") + new_header_stokes["P_int"] = (debiased_P_diluted, "Integrated polarization degree") new_header_stokes["sP_int"] = (np.ceil(P_diluted_err * 1000.0) / 1000.0, "Integrated polarization degree error") 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_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_data_mask, new_header_stokes, new_s_IQU_stat - else: - return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_data_mask, new_header_stokes + return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_Stokes_stat_cov, new_data_mask, new_header_stokes def rotate_data(data_array, error_array, data_mask, headers): diff --git a/package/src/emission_center.py b/package/src/emission_center.py index 6e3b892..a8f2092 100755 --- a/package/src/emission_center.py +++ b/package/src/emission_center.py @@ -43,7 +43,9 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None): if target is None: target = Stokes[0].header["TARGNAME"] - fig = figure(figsize=(8, 8.5), layout="constrained") + ratiox = max(int(stkI.shape[1] / (stkI.shape[0])), 1) + ratioy = max(int((stkI.shape[0]) / stkI.shape[1]), 1) + fig = figure(figsize=(8 * ratiox, 8 * ratioy), 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)