diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot.png index bad3f27..84ae9e2 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P.png index d558dc5..0c1b977 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_err.png index 4604557..5b19875 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_err.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_flux.png index ecd9edb..e8b5320 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_flux.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_flux.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRi.png index bd85ae3..cee85cf 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRi.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRp.png index f1292ab..ebd3494 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRp.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRp.png differ diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index bc2dc4f..57a7239 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -110,8 +110,8 @@ def main(): # Polarization map output figname = 'NGC1068_FOC' #target/intrument name figtype = '_combine_FWHM020_rot' #additionnal informations - SNRp_cut = 10 #P measurments with SNR>3 - SNRi_cut = 130 #I measurments with SNR>30, which implies an uncertainty in P of 4.7%. + SNRp_cut = 20. #P measurments with SNR>3 + SNRi_cut = 200 #I measurments with SNR>30, which implies an uncertainty in P of 4.7%. step_vec = 1 #plot all vectors in the array. if step_vec = 2, then every other vector will be plotted ##### Pipeline start @@ -172,7 +172,7 @@ def main(): # 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, pol_flux = proj_red.compute_Stokes(data_array, error_array, data_mask, headers, FWHM=smoothing_FWHM, scale=smoothing_scale, smoothing=smoothing_function) + I_stokes, Q_stokes, U_stokes, Stokes_cov = proj_red.compute_Stokes(data_array, error_array, data_mask, headers, FWHM=smoothing_FWHM, scale=smoothing_scale, smoothing=smoothing_function) ## Step 3: # Rotate images to have North up @@ -183,9 +183,9 @@ def main(): [np.sin(-alpha), np.cos(-alpha)]]) rectangle[0:2] = np.dot(mrot, np.asarray(rectangle[0:2]))+np.array(data_array.shape[1:])/2 rectangle[4] = alpha - I_stokes, Q_stokes, U_stokes, Stokes_cov, pol_flux, data_mask, headers = proj_red.rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, pol_flux, data_mask, headers, -ref_header['orientat'], SNRi_cut=None) + I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers = proj_red.rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, -ref_header['orientat'], 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, pol_flux, headers) + 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, headers) ## Step 4: # Save image to FITS. diff --git a/src/lib/plots.py b/src/lib/plots.py index 56798de..e5d5060 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -262,20 +262,20 @@ def polarization_map(Stokes, rectangle=None, SNRp_cut=3., SNRi_cut=30., step_vec elif display.lower() in ['pol_flux']: # Display polarisation flux pf_mask = (stkI.data > 0.) * (pol.data > 0.) - vmin, vmax = 0., np.max(stkI.data[pf_mask]*convert_flux*pol.data[pf_mask]/100.) - im = ax.imshow(stkI.data*convert_flux*pol.data/100.,extent=[-stkI.data.shape[1]/2.,stkI.data.shape[1]/2.,-stkI.data.shape[0]/2.,stkI.data.shape[0]/2.], vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + vmin, vmax = 0., np.max(stkI.data[pf_mask]*convert_flux*pol.data[pf_mask]) + im = ax.imshow(stkI.data*convert_flux*pol.data,extent=[-stkI.data.shape[1]/2.,stkI.data.shape[1]/2.,-stkI.data.shape[0]/2.,stkI.data.shape[0]/2.], vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") levelsI = np.linspace(SNRi_cut, np.max(SNRi[SNRi > 0.]), 10) cont = ax.contour(SNRi, extent=[-SNRi.shape[1]/2.,SNRi.shape[1]/2.,-SNRi.shape[0]/2.,SNRi.shape[0]/2.], levels=levelsI, colors='grey', linewidths=0.5) elif display.lower() in ['p','pol','pol_deg']: # Display polarization degree map vmin, vmax = 0., 100. - im = ax.imshow(pol.data,extent=[-pol.data.shape[1]/2.,pol.data.shape[1]/2.,-pol.data.shape[0]/2.,pol.data.shape[0]/2.], vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + im = ax.imshow(pol.data*100.,extent=[-pol.data.shape[1]/2.,pol.data.shape[1]/2.,-pol.data.shape[0]/2.,pol.data.shape[0]/2.], vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$P$ [%]") elif display.lower() in ['s_p','pol_err','pol_deg_err']: # Display polarization degree error map - vmin, vmax = 0., 5. - im = ax.imshow(pol_err.data,extent=[-pol_err.data.shape[1]/2.,pol_err.data.shape[1]/2.,-pol_err.data.shape[0]/2.,pol_err.data.shape[0]/2.], vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + vmin, vmax = 0., 10. + im = ax.imshow(pol_err.data*100.,extent=[-pol_err.data.shape[1]/2.,pol_err.data.shape[1]/2.,-pol_err.data.shape[0]/2.,pol_err.data.shape[0]/2.], vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\sigma_P$ [%]") elif display.lower() in ['snr','snri']: # Display I_stokes signal-to-noise map @@ -307,7 +307,7 @@ def polarization_map(Stokes, rectangle=None, SNRp_cut=3., SNRi_cut=30., step_vec X, Y = np.meshgrid(np.linspace(-stkI.data.shape[0]/2.,stkI.data.shape[0]/2.,stkI.data.shape[0]), np.linspace(-stkI.data.shape[1]/2.,stkI.data.shape[1]/2.,stkI.data.shape[1])) U, V = pol.data*np.cos(np.pi/2.+pang.data*np.pi/180.), pol.data*np.sin(np.pi/2.+pang.data*np.pi/180.) - Q = ax.quiver(X[::step_vec,::step_vec],Y[::step_vec,::step_vec],U[::step_vec,::step_vec],V[::step_vec,::step_vec],units='xy',angles='uv',scale=50.,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,color='w') + Q = ax.quiver(X[::step_vec,::step_vec],Y[::step_vec,::step_vec],U[::step_vec,::step_vec],V[::step_vec,::step_vec],units='xy',angles='uv',scale=0.5,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,color='w') pol_sc = AnchoredSizeBar(ax.transData, 2., r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w', fontproperties=fontprops) ax.add_artist(pol_sc) @@ -333,8 +333,8 @@ def polarization_map(Stokes, rectangle=None, SNRp_cut=3., SNRi_cut=30., step_vec IU_int_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[0,2][mask]**2)) QU_int_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[1,2][mask]**2)) - P_int = np.sqrt(Q_int**2+U_int**2)/I_int*100. - P_int_err = (100./I_int)*np.sqrt((Q_int**2*Q_int_err**2 + U_int**2*U_int_err**2 + 2.*Q_int*U_int*QU_int_err)/(Q_int**2 + U_int**2) + ((Q_int/I_int)**2 + (U_int/I_int)**2)*I_int_err**2 - 2.*(Q_int/I_int)*IQ_int_err - 2.*(U_int/I_int)*IU_int_err) + P_int = np.sqrt(Q_int**2+U_int**2)/I_int + P_int_err = (1./I_int)*np.sqrt((Q_int**2*Q_int_err**2 + U_int**2*U_int_err**2 + 2.*Q_int*U_int*QU_int_err)/(Q_int**2 + U_int**2) + ((Q_int/I_int)**2 + (U_int/I_int)**2)*I_int_err**2 - 2.*(Q_int/I_int)*IQ_int_err - 2.*(U_int/I_int)*IU_int_err) PA_int = princ_angle((90./np.pi)*np.arctan2(U_int,Q_int)) PA_int_err = (90./(np.pi*(Q_int**2 + U_int**2)))*np.sqrt(U_int**2*Q_int_err**2 + Q_int**2*U_int_err**2 - 2.*Q_int*U_int*QU_int_err) @@ -351,15 +351,15 @@ def polarization_map(Stokes, rectangle=None, SNRp_cut=3., SNRi_cut=30., step_vec IU_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[0,2]**2)) QU_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[1,2]**2)) - P_diluted = np.sqrt(Q_diluted**2+U_diluted**2)/I_diluted*100. - P_diluted_err = (100./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**2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err) + 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**2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err) #P_diluted_err = np.sqrt(2/n_pix)*100. PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted,Q_diluted)) PA_diluted_err = (90./(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.*Q_diluted*U_diluted*QU_diluted_err) #PA_diluted_err = P_diluted_err/(2.*P_diluted)*180./np.pi - ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(pivot_wav,sci_not(I_diluted*convert_flux,I_diluted_err*convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_diluted,P_diluted_err)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_diluted,PA_diluted_err), color='white', fontsize=16, xy=(0.01, 0.92), xycoords='axes fraction') + ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(pivot_wav,sci_not(I_diluted*convert_flux,I_diluted_err*convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_diluted*100.,P_diluted_err*100.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_diluted,PA_diluted_err), color='white', fontsize=16, xy=(0.01, 0.92), xycoords='axes fraction') ax.coords.grid(True, color='white', ls='dotted', alpha=0.5) ax.coords[0].set_axislabel('Right Ascension (J2000)') diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 11391cc..344390c 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -1059,9 +1059,6 @@ def compute_Stokes(data_array, error_array, data_mask, headers, +45/-45deg linear polarization intensity Stokes_cov : numpy.ndarray Covariance matrix of the Stokes parameters I, Q, U. - pol_flux : numpy.ndarray - Array containing the transmittance corrected fluxes from the multiple - polarizer plates """ # Check that all images are from the same instrument instr = headers[0]['instrume'] @@ -1115,7 +1112,7 @@ def compute_Stokes(data_array, error_array, data_mask, headers, norm = pol_eff[1]*pol_eff[2]*np.sin(-2.*theta[1]+2.*theta[2]) \ + pol_eff[2]*pol_eff[0]*np.sin(-2.*theta[2]+2.*theta[0]) \ + pol_eff[0]*pol_eff[1]*np.sin(-2.*theta[0]+2.*theta[1]) - globals()['a_stokes'] = np.zeros((3,3)) + a_stokes = np.zeros((3,3)) for i in range(3): a_stokes[0,i] = pol_eff[(i+1)%3]*pol_eff[(i+2)%3]*np.sin(-2.*theta[(i+1)%3]+2.*theta[(i+2)%3])/norm a_stokes[1,i] = (-pol_eff[(i+1)%3]*np.sin(2.*theta[(i+1)%3]) + pol_eff[(i+2)%3]*np.sin(2.*theta[(i+2)%3]))/norm @@ -1138,12 +1135,27 @@ def compute_Stokes(data_array, error_array, data_mask, headers, #Stokes covariance matrix Stokes_cov = np.zeros((3,3,I_stokes.shape[0],I_stokes.shape[1])) - Stokes_cov[0,0] = (4./9.)*(pol_cov[0,0]+pol_cov[1,1]+pol_cov[2,2]) + (8./9.)*(pol_cov[0,1]+pol_cov[0,2]+pol_cov[1,2]) - Stokes_cov[1,1] = (4./3.)*(pol_cov[1,1]+pol_cov[2,2]) - (8./3.)*pol_cov[1,2] - Stokes_cov[2,2] = (4./9.)*(4.*pol_cov[0,0]+pol_cov[1,1]+pol_cov[2,2]) - (8./3.)*(2.*pol_cov[0,1]+2.*pol_cov[0,2]-pol_cov[1,2]) - Stokes_cov[0,1] = Stokes_cov[1,0] = (4./(3.*np.sqrt(3.)))*(pol_cov[1,1]-pol_cov[2,2]+pol_cov[0,1]-pol_cov[0,2]) - Stokes_cov[0,2] = Stokes_cov[2,0] = (4./9.)*(2.*pol_cov[0,0]-pol_cov[1,1]-pol_cov[2,2]+pol_cov[0,1]+pol_cov[0,2]-2.*pol_cov[1,2]) - Stokes_cov[1,2] = Stokes_cov[2,1] = (4./(3.*np.sqrt(3.)))*(-pol_cov[1,1]+pol_cov[2,2]+2.*pol_cov[0,1]-2.*pol_cov[0,2]) + Stokes_cov[0,0] = a_stokes[0,0]**2*pol_cov[0,0]+a_stokes[0,1]**2*pol_cov[1,1]+a_stokes[0,2]**2*pol_cov[2,2] + 2*(a_stokes[0,0]*a_stokes[0,1]*pol_cov[0,1]+a_stokes[0,0]*a_stokes[0,2]*pol_cov[0,2]+a_stokes[0,1]*a_stokes[0,2]*pol_cov[1,2]) + Stokes_cov[1,1] = a_stokes[1,0]**2*pol_cov[0,0]+a_stokes[1,1]**2*pol_cov[1,1]+a_stokes[1,2]**2*pol_cov[2,2] + 2*(a_stokes[1,0]*a_stokes[1,1]*pol_cov[0,1]+a_stokes[1,0]*a_stokes[1,2]*pol_cov[0,2]+a_stokes[1,1]*a_stokes[1,2]*pol_cov[1,2]) + Stokes_cov[2,2] = a_stokes[2,0]**2*pol_cov[0,0]+a_stokes[2,1]**2*pol_cov[1,1]+a_stokes[2,2]**2*pol_cov[2,2] + 2*(a_stokes[2,0]*a_stokes[2,1]*pol_cov[0,1]+a_stokes[2,0]*a_stokes[2,2]*pol_cov[0,2]+a_stokes[2,1]*a_stokes[2,2]*pol_cov[1,2]) + Stokes_cov[0,1] = Stokes_cov[1,0] = a_stokes[0,0]*a_stokes[1,0]*pol_cov[0,0]+a_stokes[0,1]*a_stokes[1,1]*pol_cov[1,1]+a_stokes[0,2]*a_stokes[1,2]*pol_cov[2,2]+(a_stokes[0,0]*a_stokes[1,1]+a_stokes[1,0]*a_stokes[0,1])*pol_cov[0,1]+(a_stokes[0,0]*a_stokes[1,2]+a_stokes[1,0]*a_stokes[0,2])*pol_cov[0,2]+(a_stokes[0,1]*a_stokes[1,2]+a_stokes[1,1]*a_stokes[0,2])*pol_cov[1,2] + Stokes_cov[0,2] = Stokes_cov[2,0] = a_stokes[0,0]*a_stokes[2,0]*pol_cov[0,0]+a_stokes[0,1]*a_stokes[2,1]*pol_cov[1,1]+a_stokes[0,2]*a_stokes[2,2]*pol_cov[2,2]+(a_stokes[0,0]*a_stokes[2,1]+a_stokes[2,0]*a_stokes[0,1])*pol_cov[0,1]+(a_stokes[0,0]*a_stokes[2,2]+a_stokes[2,0]*a_stokes[0,2])*pol_cov[0,2]+(a_stokes[0,1]*a_stokes[2,2]+a_stokes[2,1]*a_stokes[0,2])*pol_cov[1,2] + Stokes_cov[1,2] = Stokes_cov[2,1] = a_stokes[1,0]*a_stokes[2,0]*pol_cov[0,0]+a_stokes[1,1]*a_stokes[2,1]*pol_cov[1,1]+a_stokes[1,2]*a_stokes[2,2]*pol_cov[2,2]+(a_stokes[1,0]*a_stokes[2,1]+a_stokes[2,0]*a_stokes[1,1])*pol_cov[0,1]+(a_stokes[1,0]*a_stokes[2,2]+a_stokes[2,0]*a_stokes[1,2])*pol_cov[0,2]+(a_stokes[1,1]*a_stokes[2,2]+a_stokes[2,1]*a_stokes[1,2])*pol_cov[1,2] + + C1 = 2.*pol_eff[0]*pol_eff[1]*pol_eff[2]/norm + dI_dtheta1 = C1*(np.cos(-2.*theta[2]+2.*theta[0])/pol_eff[1]*(pol_flux[1]-I_stokes) - np.cos(-2.*theta[0]+2.*theta[1])/pol_eff[2]*(pol_flux[2]-I_stokes)) + dI_dtheta2 = C1*(np.cos(-2.*theta[0]+2.*theta[1])/pol_eff[2]*(pol_flux[2]-I_stokes) - np.cos(-2.*theta[1]+2.*theta[2])/pol_eff[0]*(pol_flux[0]-I_stokes)) + dI_dtheta3 = C1*(np.cos(-2.*theta[1]+2.*theta[2])/pol_eff[0]*(pol_flux[0]-I_stokes) - np.cos(-2.*theta[2]+2.*theta[0])/pol_eff[1]*(pol_flux[1]-I_stokes)) + dQ_dtheta1 = C1*((np.cos(2.*theta[0])*pol_flux[1]-np.cos(2.*theta[0])*pol_flux[2])/(pol_eff[1]*pol_eff[2]) - (np.cos(-2.*theta[2]+2.*theta[0])/pol_eff[1]-np.cos(-2.*theta[0]+2.*theta[1])/pol_eff[2])*Q_stokes) + dQ_dtheta2 = C1*((np.cos(2.*theta[1])*pol_flux[2]-np.cos(2.*theta[1])*pol_flux[0])/(pol_eff[0]*pol_eff[2]) - (np.cos(-2.*theta[0]+2.*theta[1])/pol_eff[2]-np.cos(-2.*theta[1]+2.*theta[2])/pol_eff[0])*Q_stokes) + dQ_dtheta3 = C1*((np.cos(2.*theta[2])*pol_flux[0]-np.cos(2.*theta[2])*pol_flux[1])/(pol_eff[0]*pol_eff[1]) - (np.cos(-2.*theta[1]+2.*theta[2])/pol_eff[0]-np.cos(-2.*theta[2]+2.*theta[0])/pol_eff[1])*Q_stokes) + dU_dtheta1 = C1*((np.sin(2.*theta[0])*pol_flux[1]-np.sin(2.*theta[1])*pol_flux[2])/(pol_eff[1]*pol_eff[2]) - (np.cos(-2.*theta[2]+2.*theta[0])/pol_eff[1]-np.cos(-2.*theta[0]+2.*theta[1])/pol_eff[2])*U_stokes) + dU_dtheta2 = C1*((np.sin(2.*theta[1])*pol_flux[2]-np.sin(2.*theta[1])*pol_flux[0])/(pol_eff[0]*pol_eff[2]) - (np.cos(-2.*theta[0]+2.*theta[1])/pol_eff[2]-np.cos(-2.*theta[1]+2.*theta[2])/pol_eff[0])*U_stokes) + dU_dtheta3 = C1*((np.sin(2.*theta[2])*pol_flux[0]-np.sin(2.*theta[2])*pol_flux[1])/(pol_eff[0]*pol_eff[1]) - (np.cos(-2.*theta[1]+2.*theta[2])/pol_eff[0]-np.cos(-2.*theta[2]+2.*theta[0])/pol_eff[1])*U_stokes) + + #Stokes_cov[0,0] += (dI_dtheta1 + dI_dtheta2 + dI_dtheta3)**2*3.*np.pi/180. + #Stokes_cov[1,1] += (dQ_dtheta1 + dQ_dtheta2 + dQ_dtheta3)**2*3.*np.pi/180. + #Stokes_cov[2,2] += (dU_dtheta1 + dU_dtheta2 + dU_dtheta3)**2*3.*np.pi/180. if not(FWHM is None) and (smoothing.lower() in ['gaussian_after','gauss_after']): Stokes_array = np.array([I_stokes, Q_stokes, U_stokes]) @@ -1156,10 +1168,10 @@ def compute_Stokes(data_array, error_array, data_mask, headers, I_stokes, Q_stokes, U_stokes = Stokes_array Stokes_cov[0,0], Stokes_cov[1,1], Stokes_cov[2,2] = Stokes_error**2 - return I_stokes, Q_stokes, U_stokes, Stokes_cov, pol_flux + return I_stokes, Q_stokes, U_stokes, Stokes_cov -def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, pol_flux, headers): +def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers): """ Compute the polarization degree (in %) and angle (in deg) and their respective errors from given Stokes parameters. @@ -1176,9 +1188,6 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, pol_flux, headers): +45/-45deg linear polarization intensity Stokes_cov : numpy.ndarray Covariance matrix of the Stokes parameters I, Q, U. - pol_flux : numpy.ndarray - Array containing the transmittance corrected fluxes from the multiple - polarizer plates headers : header list List of headers corresponding to the images in data_array. ---------- @@ -1205,39 +1214,21 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, pol_flux, headers): """ #Polarization degree and angle computation I_pol = np.sqrt(Q_stokes**2 + U_stokes**2) - P = I_pol/I_stokes*100. + P = I_pol/I_stokes P[I_stokes <= 0.] = 0. PA = (90./np.pi)*np.arctan2(U_stokes,Q_stokes) - if (P>100.).any(): - print("WARNING : found pixels for which P > 100%", P[P>100.].size) + if (P>1).any(): + print("WARNING : found pixels for which P > 1", P[P>1.].size) #Associated errors - s_P = (100./I_stokes)*np.sqrt((Q_stokes**2*Stokes_cov[1,1] + U_stokes**2*Stokes_cov[2,2] + 2.*Q_stokes*U_stokes*Stokes_cov[1,2])/(Q_stokes**2 + U_stokes**2) + ((Q_stokes/I_stokes)**2 + (U_stokes/I_stokes)**2)*Stokes_cov[0,0] - 2.*(Q_stokes/I_stokes)*Stokes_cov[0,1] - 2.*(U_stokes/I_stokes)*Stokes_cov[0,2]) + s_P = (1/I_stokes)*np.sqrt((Q_stokes**2*Stokes_cov[1,1] + U_stokes**2*Stokes_cov[2,2] + 2.*Q_stokes*U_stokes*Stokes_cov[1,2])/(Q_stokes**2 + U_stokes**2) + ((Q_stokes/I_stokes)**2 + (U_stokes/I_stokes)**2)*Stokes_cov[0,0] - 2.*(Q_stokes/I_stokes)*Stokes_cov[0,1] - 2.*(U_stokes/I_stokes)*Stokes_cov[0,2]) s_PA = (90./(np.pi*(Q_stokes**2 + U_stokes**2)))*np.sqrt(U_stokes**2*Stokes_cov[1,1] + Q_stokes**2*Stokes_cov[2,2] - 2.*Q_stokes*U_stokes*Stokes_cov[1,2]) - #Error propagated from uncertainties in the direction of polarizers' axis - #uncertainty estimated to 3° (see Nota et al 1996) - k1, k2, k3 = pol_efficiency['pol0'], pol_efficiency['pol60'], pol_efficiency['pol120'] - f1, f2, f3 = pol_flux - theta1, theta2, theta3 = np.pi, np.pi/3., 2.*np.pi/3. - - norm = k2*k3*np.sin(-2.*theta2+2.*theta3) + k3*k1*np.sin(-2.*theta3+2.*theta1) + k1*k2*np.sin(-2.*theta1+2.*theta2) - C1 = 10000./(I_stokes**2*P) - C2 = P/I_stokes - dP_dtheta1 = 2.*(k1*k2*k3/norm) * (np.cos(-2.*theta3+2.*theta1)/k2 - np.cos(-2.*theta1+2.*theta2)/k3) * (((a_stokes[1,0]+a_stokes[2,0]-1.)*C1 + a_stokes[0,0]*C2)*f1 + ((a_stokes[0,1]) * (C2-C1))*f2 + ((a_stokes[0,2]) * (C2-C1))*f3) - dP_dtheta2 = 2.*(k1*k2*k3/norm) * (np.cos(-2.*theta1+2.*theta2)/k3 - np.cos(-2.*theta2+2.*theta3)/k1) * (((a_stokes[1,0]+a_stokes[2,0]-1./(1.-k3/k1*np.cos(-2.*theta2+2.*theta1)/np.cos(-2*theta1+2.*theta2)))*C1 + (a_stokes[0,0]-1./(1.-k1/k3*np.cos(-2.*theta1+2.*theta2)/np.cos(-2*theta2+2.*theta3)))*C2)*f1 + ((a_stokes[0,1]+np.cos(2.*theta2)/(a_stokes[1,2]*np.cos(2.*theta2)-a_stokes[1,1]*np.sin(2.*theta2))) * (C2-C1))*f2 + ((a_stokes[0,2]+np.sin(2.*theta2)/(a_stokes[1,2]*np.cos(2.*theta2)-a_stokes[1,1]*np.sin(2.*theta2))) * (C2-C1))*f3) - dP_dtheta3 = 2.*(k1*k2*k3/norm) * (np.cos(-2.*theta2+2.*theta3)/k1 - np.cos(-2.*theta3+2.*theta1)/k2) * (((a_stokes[1,0]+a_stokes[2,0]+1./(1.-k1/k2*np.cos(-2.*theta3+2.*theta1)/np.cos(-2*theta2+2.*theta3)))*C1 + (a_stokes[0,0]+1./(1.-k2/k1*np.cos(-2.*theta2+2.*theta3)/np.cos(-2*theta3+2.*theta1)))*C2)*f1 + ((a_stokes[0,1]+np.cos(2.*theta3)/(a_stokes[2,2]*np.cos(2.*theta3)-a_stokes[2,1]*np.sin(2.*theta3))) * (C2-C1))*f2 + ((a_stokes[0,2]+np.sin(2.*theta3)/(a_stokes[2,2]*np.cos(2.*theta3)-a_stokes[2,1]*np.sin(2.*theta3))) * (C2-C1))*f3) - - s_P_ax = np.sqrt(dP_dtheta1**2+dP_dtheta2**2+dP_dtheta3**2)*3./360. - s_PA_ax = np.ones(s_PA.shape)/np.sqrt(2)*3./360. - #Sum quadratically - s_P = np.sqrt(s_P**2 + s_P_ax**2) - s_PA = np.sqrt(s_PA**2 + s_PA_ax**2) debiased_P = np.sqrt(P**2 - s_P**2) - if (debiased_P>100.).any(): - print("WARNING : found pixels for which debiased_P > 100%", debiased_P[debiased_P>100.].size) + if (debiased_P>1.).any(): + print("WARNING : found pixels for which debiased_P > 100%", debiased_P[debiased_P>1.].size) #Compute the total exposure time so that #I_stokes*exp_tot = N_tot the total number of events @@ -1262,7 +1253,7 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, pol_flux, headers): 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, pol_flux, data_mask, headers, ang, SNRi_cut=None): +def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, ang, 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 @@ -1280,9 +1271,6 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, pol_flux, data_mask, +45/-45deg linear polarization intensity Stokes_cov : numpy.ndarray Covariance matrix of the Stokes parameters I, Q, U. - pol_flux : numpy.ndarray - Array containing the transmittance corrected fluxes from the multiple - polarizer plates headers : header list List of headers corresponding to the reduced images. ang : float @@ -1305,8 +1293,6 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, pol_flux, data_mask, accounting for +45/-45deg linear polarization intensity. new_Stokes_cov : numpy.ndarray Updated covariance matrix of the Stokes parameters I, Q, U. - new_pol_flux : numpy.ndarray - Rotated fluxes from the multiple polarizer plates new_headers : header list Updated list of headers corresponding to the reduced images accounting for the new orientation angle. @@ -1329,7 +1315,6 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, pol_flux, data_mask, new_Q_stokes = np.cos(2*alpha)*Q_stokes + np.sin(2*alpha)*U_stokes new_U_stokes = -np.sin(2*alpha)*Q_stokes + np.cos(2*alpha)*U_stokes - new_pol_flux = copy.deepcopy(pol_flux) #Compute new covariance matrix on rotated parameters new_Stokes_cov = copy.deepcopy(Stokes_cov) @@ -1345,7 +1330,6 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, pol_flux, data_mask, new_U_stokes = sc_rotate(new_U_stokes, ang, reshape=False, cval=0.) new_data_mask = sc_rotate(data_mask, ang, reshape=False, cval=True) for i in range(3): - new_pol_flux[i] = sc_rotate(new_pol_flux[i], ang, reshape=False, cval=0.) for j in range(3): new_Stokes_cov[i,j] = sc_rotate(new_Stokes_cov[i,j], ang, reshape=False, cval=0.) @@ -1381,7 +1365,6 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, pol_flux, data_mask, new_U_stokes[new_I_stokes == 0.] = 0. new_Q_stokes[np.isnan(new_Q_stokes)] = 0. new_U_stokes[np.isnan(new_U_stokes)] = 0. - new_pol_flux[np.isnan(new_pol_flux)] = 0. new_Stokes_cov[np.isnan(new_Stokes_cov)] = fmax - return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_pol_flux, new_data_mask, new_headers + return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_data_mask, new_headers