add comparison to Kishimoto's pipeline and context for M87 and IC5063

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']])))

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)