From ce927c3409246bd873246e5a73ecbd8bb3ebff9d Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Thu, 28 Mar 2024 17:31:58 +0100 Subject: [PATCH] fix integrated flux display, analysis default snri --- src/FOC_reduction.py | 4 ++-- src/analysis.py | 2 +- src/lib/plots.py | 33 ++++++++++++++++----------------- 3 files changed, 19 insertions(+), 20 deletions(-) diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 1b41ee9..21ed131 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -172,7 +172,6 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= figname = "_".join([figname, figtype]) if figtype != "" else figname Stokes_test = proj_fits.save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, headers, data_mask, figname, data_folder=data_folder, return_hdul=True) - data_mask = Stokes_test[-1].data.astype(bool) # Step 5: # crop to desired region of interest (roi) @@ -181,8 +180,9 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= stokescrop = proj_plots.crop_Stokes(deepcopy(Stokes_test), norm=LogNorm()) stokescrop.crop() stokescrop.write_to("/".join([data_folder, figname+".fits"])) - Stokes_test, data_mask, headers = stokescrop.hdul_crop, stokescrop.data_mask, [dataset.header for dataset in stokescrop.hdul_crop] + Stokes_test, headers = stokescrop.hdul_crop, [dataset.header for dataset in stokescrop.hdul_crop] + data_mask = Stokes_test['data_mask'].data.astype(bool) print("F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(headers[0]['photplam'], *sci_not( Stokes_test[0].data[data_mask].sum()*headers[0]['photflam'], np.sqrt(Stokes_test[3].data[0, 0][data_mask].sum())*headers[0]['photflam'], 2, out=int))) print("P_int = {0:.1f} ± {1:.1f} %".format(headers[0]['p_int']*100., np.ceil(headers[0]['p_int_err']*1000.)/10.)) diff --git a/src/analysis.py b/src/analysis.py index 9d2aaff..e2d2fea 100755 --- a/src/analysis.py +++ b/src/analysis.py @@ -7,7 +7,7 @@ options = "hf:p:i:l:" long_options = ["help", "fits=", "snrp=", "snri=", "lim="] fits_path = None -SNRp_cut, SNRi_cut = 3, 30 +SNRp_cut, SNRi_cut = 3, 3 flux_lim = None out_txt = None diff --git a/src/lib/plots.py b/src/lib/plots.py index 6675845..d305784 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -1408,16 +1408,16 @@ class crop_Stokes(crop_map): self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, button=[1]) # Update integrated values - mask = np.logical_and(self.hdul_crop[-1].data.astype(bool), self.hdul_crop[0].data > 0) - I_diluted = self.hdul_crop[0].data[mask].sum() - Q_diluted = self.hdul_crop[1].data[mask].sum() - U_diluted = self.hdul_crop[2].data[mask].sum() - I_diluted_err = np.sqrt(np.sum(self.hdul_crop[3].data[0, 0][mask])) - Q_diluted_err = np.sqrt(np.sum(self.hdul_crop[3].data[1, 1][mask])) - U_diluted_err = np.sqrt(np.sum(self.hdul_crop[3].data[2, 2][mask])) - IQ_diluted_err = np.sqrt(np.sum(self.hdul_crop[3].data[0, 1][mask]**2)) - IU_diluted_err = np.sqrt(np.sum(self.hdul_crop[3].data[0, 2][mask]**2)) - QU_diluted_err = np.sqrt(np.sum(self.hdul_crop[3].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)) P_diluted = np.sqrt(Q_diluted**2+U_diluted**2)/I_diluted P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted ** @@ -1436,7 +1436,7 @@ class crop_Stokes(crop_map): @property def data_mask(self): - return self.hdul_crop[-1].data.astype(int) + return self.hdul_crop['data_mask'].data.astype(int) class image_lasso_selector(object): @@ -2341,12 +2341,12 @@ class pol_map(object): if hasattr(self, 'quiver'): self.quiver.remove() self.quiver = ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=1./self.vec_scale, scale_units='xy', pivot='mid', headwidth=0., - headlength=0., headaxislength=0., width=0.15, linewidth=0.5, color='white', edgecolor='black') + headlength=0., headaxislength=0., width=0.2, linewidth=0.3, color='white', edgecolor='black') fig.canvas.draw_idle() return self.quiver else: ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=1./self.vec_scale, scale_units='xy', pivot='mid', headwidth=0., - headlength=0., headaxislength=0., width=0.15, linewidth=0.5, color='white', edgecolor='black') + headlength=0., headaxislength=0., width=0.2, linewidth=0.3, color='white', edgecolor='black') fig.canvas.draw_idle() def pol_int(self, fig=None, ax=None): @@ -2380,9 +2380,8 @@ class pol_map(object): P_cut_err = np.sqrt((Q_cut**2*Q_cut_err**2 + U_cut**2*U_cut_err**2 + 2.*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.*(Q_cut/I_cut)*IQ_cut_err - 2.*(U_cut/I_cut)*IU_cut_err)/I_cut - PA_cut = princ_angle(np.degrees((1./2.)*np.arctan2(U_cut, Q_cut))) - PA_cut_err = princ_angle(np.degrees((1./(2.*(Q_cut**2+U_cut**2)))*np.sqrt(U_cut**2 * - Q_cut_err**2 + Q_cut**2*U_cut_err**2 - 2.*Q_cut*U_cut*QU_cut_err))) + PA_cut = princ_angle((90./np.pi)*np.arctan2(U_cut, Q_cut)) + PA_cut_err = (90./(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.*Q_cut*U_cut*QU_cut_err) else: s_I = np.sqrt(self.IQU_cov[0, 0]) @@ -2424,7 +2423,7 @@ class pol_map(object): P_cut_err = np.sqrt((Q_cut**2*Q_cut_err**2 + U_cut**2*U_cut_err**2 + 2.*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.*(Q_cut/I_cut)*IQ_cut_err - 2.*(U_cut/I_cut)*IU_cut_err)/I_cut - PA_cut = 360.-princ_angle((90./np.pi)*np.arctan2(U_cut, Q_cut)) + PA_cut = princ_angle((90./np.pi)*np.arctan2(U_cut, Q_cut)) PA_cut_err = (90./(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.*Q_cut*U_cut*QU_cut_err) if hasattr(self, 'cont'):