diff --git a/plots/IC5063/IC5063_regions_pol.png b/plots/IC5063/IC5063_regions_pol.png new file mode 100644 index 0000000..4816761 Binary files /dev/null and b/plots/IC5063/IC5063_regions_pol.png differ diff --git a/plots/IC5063/Observed_polarization_spectrum_IC5063.png b/plots/IC5063/Observed_polarization_spectrum_IC5063.png new file mode 100644 index 0000000..a2cf67e Binary files /dev/null and b/plots/IC5063/Observed_polarization_spectrum_IC5063.png differ diff --git a/plots/M87/ObservationLog.png b/plots/M87/ObservationLog.png new file mode 100644 index 0000000..a20586a Binary files /dev/null and b/plots/M87/ObservationLog.png differ diff --git a/plots/NGC1068/5144/NGC1068_K_comparison.png b/plots/NGC1068/5144/NGC1068_K_comparison.png new file mode 100644 index 0000000..577a38a Binary files /dev/null and b/plots/NGC1068/5144/NGC1068_K_comparison.png differ diff --git a/plots/NGC1068/5144/NGC1068_K_pol_ang.png b/plots/NGC1068/5144/NGC1068_K_pol_ang.png new file mode 100644 index 0000000..cdbde8a Binary files /dev/null and b/plots/NGC1068/5144/NGC1068_K_pol_ang.png differ diff --git a/plots/NGC1068/5144/NGC1068_K_pol_deg.png b/plots/NGC1068/5144/NGC1068_K_pol_deg.png new file mode 100644 index 0000000..5187448 Binary files /dev/null and b/plots/NGC1068/5144/NGC1068_K_pol_deg.png differ diff --git a/plots/NGC1068/5144/NGC1068_K_pol_diff.png b/plots/NGC1068/5144/NGC1068_K_pol_diff.png new file mode 100644 index 0000000..e805c47 Binary files /dev/null and b/plots/NGC1068/5144/NGC1068_K_pol_diff.png differ diff --git a/plots/NGC1068/5144/NGC1068_K_polang_diff.png b/plots/NGC1068/5144/NGC1068_K_polang_diff.png new file mode 100644 index 0000000..8d9c67d Binary files /dev/null and b/plots/NGC1068/5144/NGC1068_K_polang_diff.png differ diff --git a/src/comparison_Kishimoto.py b/src/comparison_Kishimoto.py index dd4b6ae..ebdb05b 100755 --- a/src/comparison_Kishimoto.py +++ b/src/comparison_Kishimoto.py @@ -47,7 +47,7 @@ for d in data_K: data_K[d] = zeropad(data_K[d],shape) #shift array to get same information in same pixel -data_arr, error_ar, heads, data_msk, shifts, shifts_err = align_data(np.array([data_S['I'],data_K['I']]), [header, header], error_array=np.array([data_S['sI'],data_K['sI']]), background=np.array([bkg_S,bkg_K]), upsample_factor=10., return_shifts=True) +data_arr, error_ar, heads, data_msk, shifts, shifts_err = align_data(np.array([data_S['I'],data_K['I']]), [header, header], error_array=np.array([data_S['sI'],data_K['sI']]), background=np.array([bkg_S,bkg_K]), upsample_factor=10., ref_center='center', return_shifts=True) for d in data_K: data_K[d] = shift(data_K[d],shifts[1],order=1,cval=0.) @@ -57,9 +57,10 @@ for d in [data_S, data_K]: d[i][np.isnan(d[i])] = 0. d['P'] = np.where(np.logical_and(np.isfinite(d['I']),d['I']>0.),np.sqrt(d['Q']**2+d['U']**2)/d['I'],0.) d['sP'] = np.where(np.logical_and(np.isfinite(d['I']),d['I']>0.),np.sqrt((d['Q']**2*d['sQ']**2+d['U']**2*d['sU']**2)/(d['Q']**2+d['U']**2)+((d['Q']/d['I'])**2+(d['U']/d['I'])**2)*d['sI']**2)/d['I'],0.) + d['d_P'] = np.where(np.logical_and(np.isfinite(d['P']),np.isfinite(d['sP'])),np.sqrt(d['P']**2-d['sP']**2),0.) d['PA'] = 0.5*np.arctan2(d['U'],d['Q'])+np.pi - d['SNRp'] = np.zeros(d['P'].shape) - d['SNRp'][d['sP']>0.] = d['P'][d['sP']>0.]/d['sP'][d['sP']>0.] + d['SNRp'] = np.zeros(d['d_P'].shape) + d['SNRp'][d['sP']>0.] = d['d_P'][d['sP']>0.]/d['sP'][d['sP']>0.] d['SNRi'] = np.zeros(d['I'].shape) d['SNRi'][d['sI']>0.] = d['I'][d['sI']>0.]/d['sI'][d['sI']>0.] d['mask'] = np.logical_and(d['SNRi']>SNRi_cut,d['SNRp']>SNRp_cut) @@ -73,7 +74,7 @@ bins=int(data_S['mask'].sum()/5) bin_size=1./bins mod_p = np.linspace(0.,1.,300) for d in [data_S, data_K]: - d['hist'], d['bin_edges'] = np.histogram(d['P'][d['mask']],bins=bins,range=(0.,1.)) + d['hist'], d['bin_edges'] = np.histogram(d['d_P'][d['mask']],bins=bins,range=(0.,1.)) d['binning'] = bin_centers(d['bin_edges']) peak, bins_fwhm = d['binning'][np.argmax(d['hist'])], d['binning'][d['hist']>d['hist'].max()/2.] fwhm = bins_fwhm[1]-bins_fwhm[0] @@ -92,7 +93,7 @@ ax_p.errorbar(data_K['binning'],data_K['hist'],xerr=bin_size/2.,fmt='r.',ecolor= ax_p.plot(mod_p,gauss(mod_p,*data_K['hist_popt']),'r--',label='mean = {1:.2f}, stdev = {2:.2f}'.format(*data_K['hist_popt'])) ax_p.set(xlabel="Polarization degree",ylabel="Counts",title="Histogram of polarization degree computed in the cut for both pipelines.") ax_p.legend() -fig_p.savefig(path_join(root_dir_plot_S,"NGC1068_K_pol_deg.png"),bbox_inches="tight") +fig_p.savefig(path_join(root_dir_plot_S,"NGC1068_K_pol_deg.png"),bbox_inches="tight",dpi=300) ##### ###Compute angular difference between the maps in cut @@ -106,7 +107,35 @@ im_pa = ax_pa.imshow(np.cos(2*dtheta), vmin=-1., vmax=1., origin='lower', cmap=' cbar_pa = plt.colorbar(im_pa, cax=cbar_ax_pa, label=r"$\zeta = \cos\left( 2 \cdot \delta\theta_P \right)$") ax_pa.coords[0].set_axislabel('Right Ascension (J2000)') ax_pa.coords[1].set_axislabel('Declination (J2000)') -fig_pa.savefig(path_join(root_dir_plot_S,"NGC1068_K_pol_ang.png"),bbox_inches="tight") +fig_pa.savefig(path_join(root_dir_plot_S,"NGC1068_K_pol_ang.png"),bbox_inches="tight",dpi=300) + +##### +###Compute power uncertainty difference between the maps in cut +##### +eta = np.where(data_S['mask'], np.abs(data_K['d_P']-data_S['d_P'])/np.sqrt(data_S['sP']**2+data_K['sP']**2)/2., np.nan) +fig_dif_p = plt.figure(num="Polarization power difference ratio") +ax_dif_p = fig_dif_p.add_subplot(111, projection=wcs) +cbar_ax_dif_p = fig_dif_p.add_axes([0.88, 0.12, 0.01, 0.75]) +ax_dif_p.set_title(r"Degree of difference $\eta$ of the polarization from the 2 pipelines in the cut") +im_dif_p = ax_dif_p.imshow(eta, vmin=0., vmax=2., origin='lower', cmap='bwr_r', label=r"$\eta$ between this pipeline and Kishimoto's") +cbar_dif_p = plt.colorbar(im_dif_p, cax=cbar_ax_dif_p, label=r"$\eta = \frac{2 \left|P^K-P^S\right|}{\sqrt{{\sigma^K_P}^2+{\sigma^S_P}^2}}$") +ax_dif_p.coords[0].set_axislabel('Right Ascension (J2000)') +ax_dif_p.coords[1].set_axislabel('Declination (J2000)') +fig_dif_p.savefig(path_join(root_dir_plot_S,"NGC1068_K_pol_diff.png"),bbox_inches="tight",dpi=300) + +##### +###Compute angle uncertainty difference between the maps in cut +##### +eta = np.where(data_S['mask'], np.abs(data_K['PA']-data_S['PA'])/np.sqrt(data_S['sPA']**2+data_K['sPA']**2)/2., np.nan) +fig_dif_pa = plt.figure(num="Polarization angle difference ratio") +ax_dif_pa = fig_dif_pa.add_subplot(111, projection=wcs) +cbar_ax_dif_pa = fig_dif_pa.add_axes([0.88, 0.12, 0.01, 0.75]) +ax_dif_pa.set_title(r"Degree of difference $\eta$ of the polarization from the 2 pipelines in the cut") +im_dif_pa = ax_dif_pa.imshow(eta, vmin=0., vmax=2., origin='lower', cmap='bwr_r', label=r"$\eta$ between this pipeline and Kishimoto's") +cbar_dif_pa = plt.colorbar(im_dif_pa, cax=cbar_ax_dif_pa, label=r"$\eta = \frac{2 \left|\theta_P^K-\theta_P^S\right|}{\sqrt{{\sigma^K_{\theta_P}}^2+{\sigma^S_{\theta_P}}^2}}$") +ax_dif_pa.coords[0].set_axislabel('Right Ascension (J2000)') +ax_dif_pa.coords[1].set_axislabel('Declination (J2000)') +fig_dif_pa.savefig(path_join(root_dir_plot_S,"NGC1068_K_polang_diff.png"),bbox_inches="tight",dpi=300) ##### ###display both polarization maps to check consistency @@ -118,7 +147,7 @@ cbar_ax = fig.add_axes([0.88, 0.12, 0.01, 0.75]) for d in [data_S, data_K]: d['X'], d['Y'] = np.meshgrid(np.arange(d['I'].shape[1]), np.arange(d['I'].shape[0])) - d['xy_U'], d['xy_V'] = np.where(d['mask'],d['P']*np.cos(np.pi/2.+d['PA']), np.nan), np.where(d['mask'],d['P']*np.sin(np.pi/2.+d['PA']), np.nan) + d['xy_U'], d['xy_V'] = np.where(d['mask'],d['d_P']*np.cos(np.pi/2.+d['PA']), np.nan), np.where(d['mask'],d['d_P']*np.sin(np.pi/2.+d['PA']), np.nan) im0 = ax.imshow(data_S['I']*convert_flux,norm=LogNorm(data_S['I'][data_S['I']>0].min()*convert_flux,data_S['I'][data_S['I']>0].max()*convert_flux),origin='lower',cmap='gray',label=r"$I_{STOKES}$ through this pipeline") quiv0 = ax.quiver(data_S['X'],data_S['Y'],data_S['xy_U'],data_S['xy_V'],units='xy',angles='uv',scale=0.5,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.2,color='b',alpha=0.75, label="PA through this pipeline") @@ -137,7 +166,7 @@ ax.coords[1].set_ticklabel_position('l') cbar = plt.colorbar(im0, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") #plt.rcParams.update({'font.size': 8}) ax.legend(loc='upper right') -fig.savefig(path_join(root_dir_plot_S,"NGC1068_K_comparison.png"),bbox_inches="tight") +fig.savefig(path_join(root_dir_plot_S,"NGC1068_K_comparison.png"),bbox_inches="tight",dpi=300) #compute integrated polarization parameters on a specific cut for d in [data_S, data_K]: @@ -150,10 +179,11 @@ for d in [data_S, data_K]: d['P_dil'] = np.sqrt(d['Q_dil']**2+d['U_dil']**2)/d['I_dil'] d['sP_dil'] = np.sqrt((d['Q_dil']**2*d['sQ_dil']**2+d['U_dil']**2*d['sU_dil']**2)/(d['Q_dil']**2+d['U_dil']**2)+((d['Q_dil']/d['I_dil'])**2+(d['U_dil']/d['I_dil'])**2)*d['sI_dil']**2)/d['I_dil'] + d['d_P_dil'] = np.sqrt(d['P_dil']**2-d['sP_dil']**2) d['PA_dil'] = princ_angle((90./np.pi)*np.arctan2(d['U_dil'],d['Q_dil'])) d['sPA_dil'] = princ_angle((90./(np.pi*(d['Q_dil']**2+d['U_dil']**2)))*np.sqrt(d['Q_dil']**2*d['sU_dil']**2+d['U_dil']**2*d['sU_dil']**2)) -print('From this pipeline :\n', "P = {0:.2f} ± {1:.2f} %\n".format(data_S['P_dil']*100.,data_S['sP_dil']*100.), "PA = {0:.2f} ± {1:.2f} °".format(data_S['PA_dil'],data_S['sPA_dil'])) -print("From Kishimoto's pipeline :\n", "P = {0:.2f} ± {1:.2f} %\n".format(data_K['P_dil']*100.,data_K['sP_dil']*100.), "PA = {0:.2f} ± {1:.2f} °".format(data_K['PA_dil'],data_K['sPA_dil'])) +print('From this pipeline :\n', "P = {0:.2f} ± {1:.2f} %\n".format(data_S['d_P_dil']*100.,data_S['sP_dil']*100.), "PA = {0:.2f} ± {1:.2f} °".format(data_S['PA_dil'],data_S['sPA_dil'])) +print("From Kishimoto's pipeline :\n", "P = {0:.2f} ± {1:.2f} %\n".format(data_K['d_P_dil']*100.,data_K['sP_dil']*100.), "PA = {0:.2f} ± {1:.2f} °".format(data_K['PA_dil'],data_K['sPA_dil'])) #compare different types of error print("This pipeline : average sI/I={0:.2f} ; sQ/Q={1:.2f} ; sU/U={2:.2f} ; sP/P={3:.2f}".format(np.mean(data_S['sI'][data_S['mask']]/data_S['I'][data_S['mask']]),np.mean(data_S['sQ'][data_S['mask']]/data_S['Q'][data_S['mask']]),np.mean(data_S['sU'][data_S['mask']]/data_S['U'][data_S['mask']]),np.mean(data_S['sP'][data_S['mask']]/data_S['P'][data_S['mask']]))) diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 5e87c12..aacfc87 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -68,6 +68,9 @@ globals()['pol_efficiency'] = {'pol0' : 0.92, 'pol60' : 0.92, 'pol120' : 0.91} globals()['theta'] = np.array([180.*np.pi/180., 60.*np.pi/180., 120.*np.pi/180.]) # Uncertainties on the orientation of the polarizers' axes taken to be 3deg (see Nota et. al 1996, p36; Robinson & Thomson 1995) globals()['sigma_theta'] = np.array([3.*np.pi/180., 3.*np.pi/180., 3.*np.pi/180.]) +# Image shift between polarizers as measured by Hodge (1995) +globals()['pol_shift'] = {'pol0' : np.array([0.,0.])*1., 'pol60' : np.array([3.63,-0.68])*1., 'pol120' : np.array([0.65,0.20])*1.} +globals()['sigma_shift'] = {'pol0' : [0.3,0.3], 'pol60' : [0.3,0.3], 'pol120' : [0.3,0.3]} def princ_angle(ang): @@ -598,7 +601,6 @@ def rebin_array(data_array, error_array, headers, pxsize, scale, Dxy = image.shape/new_shape if (Dxy < 1.).any(): raise ValueError("Requested pixel size is below resolution.") - new_shape = np.ceil(image.shape/Dxy).astype(int) # Rebin data rebin_data = bin_ndarray(image, new_shape=new_shape, @@ -771,6 +773,10 @@ def align_data(data_array, headers, error_array=None, background=None, if do_shift: rescaled_image[i] = sc_shift(rescaled_image[i], shift, order=1, cval=0.) rescaled_error[i] = sc_shift(rescaled_error[i], shift, order=1, cval=background[i]) + else: + shift = pol_shift[headers[i]['filtnam1'].lower()] + rescaled_image[i] = sc_shift(rescaled_image[i], shift, order=1, cval=0.) + rescaled_error[i] = sc_shift(rescaled_error[i], shift, order=1, cval=background[i]) curr_mask = sc_shift(res_mask, shift, order=1, cval=False) mask_vertex = clean_ROI(curr_mask)