diff --git a/data/IC5063_x3nl030/IR/Figure_3.png b/data/IC5063_x3nl030/IR/Figure_3.png deleted file mode 100644 index a013331..0000000 Binary files a/data/IC5063_x3nl030/IR/Figure_3.png and /dev/null differ diff --git a/data/IC5063_x3nl030/IR/IC5063_WFPC2_F606W.png b/data/IC5063_x3nl030/IR/IC5063_WFPC2_F606W.png new file mode 100644 index 0000000..378c102 Binary files /dev/null and b/data/IC5063_x3nl030/IR/IC5063_WFPC2_F606W.png differ diff --git a/plots/IC5063_x3nl030/18GHz_overplot_forced.png b/plots/IC5063_x3nl030/18GHz_overplot_forced.png index 4e0bbc7..0cbc494 100644 Binary files a/plots/IC5063_x3nl030/18GHz_overplot_forced.png and b/plots/IC5063_x3nl030/18GHz_overplot_forced.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png index 01c985d..8d24221 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png index 3b3ddca..6f55929 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png index baaf05b..ad7918c 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png index a224f4d..1f08e5b 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png index d0d94ba..b83a0aa 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png index dfdb6c4..81bda65 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png index 0ccf1ac..a20f7b7 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png index e383845..79e6808 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol.png new file mode 100644 index 0000000..b943d77 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_IQU.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_IQU.png new file mode 100644 index 0000000..6f55929 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_IQU.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_I_err.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_I_err.png new file mode 100644 index 0000000..7046d45 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_I_err.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_P.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_P.png new file mode 100644 index 0000000..12c38f8 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_P.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_P_err.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_P_err.png new file mode 100644 index 0000000..4364d72 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_P_err.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_P_flux.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_P_flux.png new file mode 100644 index 0000000..bad3976 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_P_flux.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_SNRi.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_SNRi.png new file mode 100644 index 0000000..caa650b Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_SNRi.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_SNRp.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_SNRp.png new file mode 100644 index 0000000..7eea0e6 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol_SNRp.png differ diff --git a/plots/IC5063_x3nl030/IR_overplot_forced.png b/plots/IC5063_x3nl030/IR_overplot_forced.png index 5768274..6ff719a 100644 Binary files a/plots/IC5063_x3nl030/IR_overplot_forced.png and b/plots/IC5063_x3nl030/IR_overplot_forced.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010.png new file mode 100644 index 0000000..df18c0c Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_IQU.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_IQU.png new file mode 100644 index 0000000..b652894 Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_IQU.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_I_err.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_I_err.png new file mode 100644 index 0000000..ecadec8 Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_I_err.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_P.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_P.png new file mode 100644 index 0000000..fc071c4 Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_P.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_P_err.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_P_err.png new file mode 100644 index 0000000..5cb8296 Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_P_err.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_P_flux.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_P_flux.png new file mode 100644 index 0000000..f6cb02c Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_P_flux.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_SNRi.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_SNRi.png new file mode 100644 index 0000000..11458aa Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_SNRi.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_SNRp.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_SNRp.png new file mode 100644 index 0000000..243e2d3 Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM010_SNRp.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020.png new file mode 100644 index 0000000..529c0b6 Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_IQU.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_IQU.png new file mode 100644 index 0000000..b271f60 Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_IQU.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_I_err.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_I_err.png new file mode 100644 index 0000000..b24b20a Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_I_err.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_P.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_P.png new file mode 100644 index 0000000..2f6f541 Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_P.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_P_err.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_P_err.png new file mode 100644 index 0000000..11c1bf3 Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_P_err.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_P_flux.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_P_flux.png new file mode 100644 index 0000000..a581696 Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_P_flux.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_SNRi.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_SNRi.png new file mode 100644 index 0000000..e08604a Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_SNRi.png differ diff --git a/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_SNRp.png b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_SNRp.png new file mode 100644 index 0000000..6902ba1 Binary files /dev/null and b/plots/MKN463_x2rp030/MRK463_FOC_combine_FWHM020_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_Capetti_vs_us.png b/plots/NGC1068_x274020/NGC1068_Capetti_vs_us.png new file mode 100644 index 0000000..6c39f68 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_Capetti_vs_us.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_analysis.png b/plots/NGC1068_x274020/NGC1068_FOC_analysis.png new file mode 100644 index 0000000..3939435 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_analysis.png differ diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 8ca2624..523c2f7 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -20,17 +20,17 @@ from astropy.wcs import WCS def main(): ##### User inputs ## Input and output locations - globals()['data_folder'] = "../data/NGC1068_x274020/" - infiles = ['x274020at.c0f.fits','x274020bt.c0f.fits','x274020ct.c0f.fits', - 'x274020dt.c0f.fits','x274020et.c0f.fits','x274020ft.c0f.fits', - 'x274020gt.c0f.fits','x274020ht.c0f.fits','x274020it.c0f.fits'] - psf_file = 'NGC1068_f253m00.fits' - globals()['plots_folder'] = "../plots/NGC1068_x274020/" +# globals()['data_folder'] = "../data/NGC1068_x274020/" +# infiles = ['x274020at.c0f.fits','x274020bt.c0f.fits','x274020ct.c0f.fits', +# 'x274020dt.c0f.fits','x274020et.c0f.fits','x274020ft.c0f.fits', +# 'x274020gt.c0f.fits','x274020ht.c0f.fits','x274020it.c0f.fits'] +## psf_file = 'NGC1068_f253m00.fits' +# globals()['plots_folder'] = "../plots/NGC1068_x274020/" -# globals()['data_folder'] = "../data/IC5063_x3nl030/" -# infiles = ['x3nl0301r_c0f.fits','x3nl0302r_c0f.fits','x3nl0303r_c0f.fits'] + globals()['data_folder'] = "../data/IC5063_x3nl030/" + infiles = ['x3nl0301r_c0f.fits','x3nl0302r_c0f.fits','x3nl0303r_c0f.fits'] # psf_file = 'IC5063_f502m00.fits' -# globals()['plots_folder'] = "../plots/IC5063_x3nl030/" + globals()['plots_folder'] = "../plots/IC5063_x3nl030/" # globals()['data_folder'] = "../data/NGC1068_x14w010/" # infiles = ['x14w0101t_c0f.fits','x14w0102t_c0f.fits','x14w0103t_c0f.fits', @@ -122,10 +122,10 @@ def main(): crop = False #Crop to desired ROI final_display = False # Polarization map output - figname = 'NGC1068_FOC' #target/intrument name + figname = 'IC5063_FOC' #target/intrument name figtype = '_combine_FWHM020' #additionnal informations - SNRp_cut = 5. #P measurments with SNR>3 - SNRi_cut = 50. #I measurments with SNR>30, which implies an uncertainty in P of 4.7%. + SNRp_cut = 3. #P measurments with SNR>3 + SNRi_cut = 60. #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 # if step_vec = 0 then all vectors are displayed at full length @@ -176,13 +176,13 @@ def main(): # 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 = 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 if rotate_stokes: alpha = headers[0]['orientat'] I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, data_mask = proj_red.rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, -alpha, 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, headers) diff --git a/src/lib/plots.py b/src/lib/plots.py index f0e465c..70cbc9f 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -271,7 +271,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c if display is None: # If no display selected, show intensity map vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux) - im = ax.imshow(stkI.data*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + im = ax.imshow(stkI.data*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") levelsI = np.linspace(vmax*0.01, vmax*0.99, 10) print("Total flux contour levels : ", levelsI) @@ -281,7 +281,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c # 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]) - im = ax.imshow(stkI.data*convert_flux*pol.data, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + im = ax.imshow(stkI.data*convert_flux*pol.data, vmin=vmin, vmax=vmax, aspect='equal', 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}$]") levelsPf = np.linspace(vmax*0.01, vmax*0.99, 10) print("Polarized flux contour levels : ", levelsPf) @@ -290,24 +290,24 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c elif display.lower() in ['p','pol','pol_deg']: # Display polarization degree map vmin, vmax = 0., 100. - im = ax.imshow(pol.data*100., vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + im = ax.imshow(pol.data*100., vmin=vmin, vmax=vmax, aspect='equal', 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., 10. p_err = deepcopy(pol_err.data) p_err[p_err > vmax/100.] = np.nan - im = ax.imshow(p_err*100., vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + im = ax.imshow(p_err*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\sigma_P$ [%]") elif display.lower() in ['s_i','i_err']: # Display intensity error map vmin, vmax = 0., np.max(np.sqrt(stk_cov.data[0,0][stk_cov.data[0,0] > 0.])*convert_flux) - im = ax.imshow(np.sqrt(stk_cov.data[0,0])*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + im = ax.imshow(np.sqrt(stk_cov.data[0,0])*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") elif display.lower() in ['snr','snri']: # Display I_stokes signal-to-noise map vmin, vmax = 0., np.max(SNRi[SNRi > 0.]) - im = ax.imshow(SNRi, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + im = ax.imshow(SNRi, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$I_{Stokes}/\sigma_{I}$") levelsSNRi = np.linspace(SNRi_cut, vmax*0.99, 10) print("SNRi contour levels : ", levelsSNRi) @@ -316,7 +316,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c elif display.lower() in ['snrp']: # Display polarization degree signal-to-noise map vmin, vmax = SNRp_cut, np.max(SNRp[SNRp > 0.]) - im = ax.imshow(SNRp, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + im = ax.imshow(SNRp, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$P/\sigma_{P}$") levelsSNRp = np.linspace(SNRp_cut, vmax*0.99, 10) print("SNRp contour levels : ", levelsSNRp) @@ -325,7 +325,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c else: # Defaults to intensity map vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux*2.) - im = ax.imshow(stkI.data*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + im = ax.imshow(stkI.data*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") if (display is None) or not(display.lower() in ['default']): @@ -398,7 +398,7 @@ class align_maps(object): elif self.wcs_map.naxis == 3: self.wcs_map = WCS(self.map[0],naxis=[1,2]).deepcopy() self.map[0].data = self.map[0].data[1] - + self.wcs_other = WCS(self.other_map[0]).deepcopy() if self.wcs_other.naxis == 4: self.wcs_other = WCS(self.other_map[0],naxis=[1,2]).deepcopy() @@ -406,7 +406,7 @@ class align_maps(object): elif self.wcs_other.naxis == 3: self.wcs_other = WCS(self.other_map[0],naxis=[1,2]).deepcopy() self.other_map[0].data = self.other_map[0].data[1] - + try: convert_flux = self.map[0].header['photflam'] except KeyError: @@ -415,7 +415,7 @@ class align_maps(object): other_convert = self.other_map[0].header['photflam'] except KeyError: other_convert = 1. - + #Get data data = self.map[0].data other_data = self.other_map[0].data @@ -427,13 +427,13 @@ class align_maps(object): self.ax1.set_facecolor('k') vmin, vmax = 0., np.max(data[data > 0.]*convert_flux) - im1 = self.ax1.imshow(data*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + im1 = self.ax1.imshow(data*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) fontprops = fm.FontProperties(size=16) px_size = self.wcs_map.wcs.get_cdelt()[0]*3600. px_sc = AnchoredSizeBar(self.ax1.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w', fontproperties=fontprops) self.ax1.add_artist(px_sc) - + try: north_dir1 = AnchoredDirectionArrows(self.ax1.transAxes, "E", "N", length=-0.08, fontsize=0.025, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, back_length=0., head_length=10., head_width=10., angle=-self.map[0].header['orientat'], color='white', text_props={'ec': None, 'fc': 'w', 'alpha': 1, 'lw': 0.4}, arrow_props={'ec': None,'fc':'w','alpha': 1,'lw': 1}) self.ax1.add_artist(north_dir1) @@ -454,14 +454,14 @@ class align_maps(object): test = kwargs[key] except KeyError: for key_i, val_i in value: - kwargs[key_i] = val_i - im2 = self.ax2.imshow(other_data*other_convert, aspect='auto', **kwargs) + kwargs[key_i] = val_i + im2 = self.ax2.imshow(other_data*other_convert, aspect='equal', **kwargs) fontprops = fm.FontProperties(size=16) px_size = self.wcs_other.wcs.get_cdelt()[0]*3600. px_sc = AnchoredSizeBar(self.ax2.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w', fontproperties=fontprops) self.ax2.add_artist(px_sc) - + try: north_dir2 = AnchoredDirectionArrows(self.ax2.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.other_map[0].header['orientat'], color='w', arrow_props={'ec': None, 'fc': 'w', 'alpha': 1,'lw': 2}) self.ax2.add_artist(north_dir2) @@ -477,7 +477,7 @@ class align_maps(object): self.bapply = Button(self.axapply, 'Apply reference') self.axreset = self.fig.add_axes([0.60, 0.01, 0.1, 0.04]) self.breset = Button(self.axreset, 'Leave as is') - + def get_aligned_wcs(self): return self.wcs_map, self.wcs_other @@ -489,14 +489,14 @@ class align_maps(object): self.cr_map.set(data=[x,y]) self.fig.canvas.draw_idle() - + if (event.inaxes is not None) and (event.inaxes == self.ax2): x = event.xdata y = event.ydata self.cr_other.set(data=[x,y]) self.fig.canvas.draw_idle() - + def reset_align(self, event): self.wcs_map.wcs.crpix = WCS(self.map[0].header).wcs.crpix[:2] self.wcs_other.wcs.crpix = WCS(self.other_map[0].header).wcs.crpix[:2] @@ -504,7 +504,7 @@ class align_maps(object): if self.aligned: plt.close() - + self.aligned = True def apply_align(self, event): @@ -515,13 +515,13 @@ class align_maps(object): if self.aligned: plt.close() - + self.aligned = True - + def on_close_align(self, event): self.aligned = True #print(self.get_aligned_wcs()) - + def align(self): self.fig.canvas.draw() self.fig.canvas.mpl_connect('button_press_event', self.onclick_ref) @@ -547,7 +547,7 @@ class overplot_radio(align_maps): pol = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='Pol_deg_debiased' for i in range(len(self.Stokes_UV))])] pol_err = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='Pol_deg_err' for i in range(len(self.Stokes_UV))])] pang = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='Pol_ang' for i in range(len(self.Stokes_UV))])] - + other_data = self.other_map[0].data other_convert = 1. other_unit = self.other_map[0].header['bunit'] @@ -555,7 +555,7 @@ class overplot_radio(align_maps): other_unit = r"mJy/Beam" other_convert = 1e3 other_freq = self.other_map[0].header['crval3'] - + convert_flux = self.Stokes_UV[0].header['photflam'] #Compute SNR and apply cuts @@ -575,7 +575,7 @@ class overplot_radio(align_maps): #Display UV intensity map with polarization vectors vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux) - im = self.ax.imshow(stkI.data*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) + im = self.ax.imshow(stkI.data*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) cbar_ax = self.fig2.add_axes([0.95, 0.12, 0.01, 0.75]) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") @@ -600,12 +600,12 @@ class overplot_radio(align_maps): north_dir = AnchoredDirectionArrows(self.ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color='w', arrow_props={'ec': None, 'fc': 'w', 'alpha': 1,'lw': 2}) self.ax.add_artist(north_dir) - + if not(savename is None): self.fig2.savefig(savename,bbox_inches='tight',dpi=200) self.fig2.canvas.draw() - + def plot(self, levels, SNRp_cut=3., SNRi_cut=30., savename=None) -> None: self.align() if self.aligned: @@ -628,9 +628,9 @@ class overplot_pol(align_maps): pol = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='Pol_deg_debiased' for i in range(len(self.Stokes_UV))])] pol_err = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='Pol_deg_err' for i in range(len(self.Stokes_UV))])] pang = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='Pol_ang' for i in range(len(self.Stokes_UV))])] - + convert_flux = self.Stokes_UV[0].header['photflam'] - + other_data = self.other_map[0].data try: other_convert = self.other_map[0].header['photflam'] @@ -654,9 +654,9 @@ class overplot_pol(align_maps): #Display Stokes I as contours levels_stkI = np.rint(np.linspace(10,99,10))/100.*np.max(stkI.data[stkI.data > 0.]*convert_flux) - cont_stkI = self.ax.contour(stkI.data*convert_flux, transform=self.ax.get_transform(self.wcs_UV), levels=levels_stkI, colors='grey') - self.ax.clabel(cont_stkI, inline=True, fontsize=8) - + cont_stkI = self.ax.contour(stkI.data*convert_flux, transform=self.ax.get_transform(self.wcs_UV), levels=levels_stkI, colors='grey', alpha=0.) + #self.ax.clabel(cont_stkI, inline=True, fontsize=8) + self.ax.autoscale(False) #Display full size polarization vectors @@ -673,27 +673,27 @@ class overplot_pol(align_maps): test = kwargs[key] except KeyError: for key_i, val_i in value: - kwargs[key_i] = val_i + kwargs[key_i] = val_i im = self.ax.imshow(other_data*other_convert, transform=self.ax.get_transform(self.wcs_other), alpha=1., **kwargs) cbar_ax = self.fig2.add_axes([0.95, 0.12, 0.01, 0.75]) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") #Display pixel scale and North direction fontprops = fm.FontProperties(size=16) - px_size = self.wcs_other.wcs.get_cdelt()[0]*3600. + px_size = self.wcs_UV.wcs.get_cdelt()[0]*3600. px_sc = AnchoredSizeBar(self.ax.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w', fontproperties=fontprops) self.ax.add_artist(px_sc) north_dir = AnchoredDirectionArrows(self.ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color='w', arrow_props={'ec': None, 'fc': 'w', 'alpha': 1,'lw': 2}) self.ax.add_artist(north_dir) - + self.ax.set(xlabel="Right Ascension (J2000)", ylabel="Declination (J2000)", title="{0:s} overplotted with polarization vectors and Stokes I contours from HST/FOC".format(obj)) if not(savename is None): self.fig2.savefig(savename,bbox_inches='tight',dpi=200) self.fig2.canvas.draw() - + def plot(self, SNRp_cut=3., SNRi_cut=30., savename=None, **kwargs) -> None: self.align() if self.aligned: @@ -711,7 +711,7 @@ class crop_map(object): self.hdul = hdul self.header = deepcopy(self.hdul[0].header) self.wcs = WCS(self.header).deepcopy() - + self.data = deepcopy(self.hdul[0].data) try: self.convert_flux = self.header['photflam'] @@ -758,16 +758,16 @@ class crop_map(object): wcs = self.wcs if convert_flux is None: convert_flux = self.convert_flux - + vmin, vmax = 0., np.max(data[data > 0.]*convert_flux) if hasattr(self, 'im'): self.im.remove() - self.im = self.ax.imshow(data*convert_flux, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=self.mask_alpha, origin='lower') + self.im = self.ax.imshow(data*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=self.mask_alpha, origin='lower') if hasattr(self, 'cr'): self.cr[0].set_data(*wcs.wcs.crpix) else: self.cr = self.ax.plot(*wcs.wcs.crpix, 'r+') - self.fig.canvas.draw_idle() + self.fig.canvas.draw_idle() return self.im @property @@ -811,10 +811,10 @@ class crop_map(object): header = self.header data = self.data wcs = self.wcs - + vertex = self.RSextent.astype(int) shape = vertex[1::2] - vertex[0::2] - + extent = np.array(self.im.get_extent()) shape_im = extent[1::2] - extent[0::2] if (shape_im.astype(int) != shape).any() and (self.RSextent != self.extent).any(): @@ -827,7 +827,7 @@ class crop_map(object): else: self.wcs_crop.wcs.crval = wcs.wcs_pix2world([self.RScenter],1)[0] self.wcs_crop.wcs.crpix = self.RScenter-self.RSextent[::2] - + # Crop dataset self.data_crop = deepcopy(data[vertex[2]:vertex[3], vertex[0]:vertex[1]]) @@ -853,7 +853,7 @@ class crop_map(object): self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, drawtype='box', button=[1], interactive=True) self.fig.canvas.draw_idle() - + def on_close(self, event) -> None: if not hasattr(self, 'hdul_crop'): self.hdul_crop = self.hdul @@ -890,7 +890,7 @@ class crop_Stokes(crop_map): hdul = self.hdul data = self.data wcs = self.wcs - + vertex = self.RSextent.astype(int) shape = vertex[1::2] - vertex[0::2] @@ -907,7 +907,7 @@ class crop_Stokes(crop_map): else: self.wcs_crop.wcs.crval = wcs.wcs_pix2world([self.RScenter],1)[0] self.wcs_crop.wcs.crpix = self.RScenter-self.RSextent[::2] - + # Crop dataset for dataset in self.hdul_crop: if dataset.header['datatype']=='IQU_cov_matrix': @@ -941,7 +941,7 @@ class crop_Stokes(crop_map): self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, drawtype='box', button=[1], interactive=True) self.fig.canvas.draw_idle() - + @property def data_mask(self): return self.hdul_crop[-1].data @@ -967,13 +967,13 @@ class image_lasso_selector: self.ax = ax self.mask_alpha = 0.1 self.embedded = True - self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect='auto', cmap='inferno',alpha=self.mask_alpha) + self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect='equal', cmap='inferno',alpha=self.mask_alpha) plt.ion() - + lineprops = {'color': 'grey', 'linewidth': 1, 'alpha': 0.8} self.lasso = LassoSelector(self.ax, self.onselect,props=lineprops, useblit=False) self.lasso.set_visible(True) - + pix_x = np.arange(self.img.shape[0]) pix_y = np.arange(self.img.shape[1]) xv, yv = np.meshgrid(pix_y,pix_x) @@ -993,10 +993,10 @@ class image_lasso_selector: p = Path(verts) self.indices = p.contains_points(self.pix, radius=0).reshape(self.img.shape[:2]) self.update_mask() - + def update_mask(self): self.displayed.remove() - self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect='auto', cmap='inferno',alpha=self.mask_alpha) + self.displayed = self.ax.imshow(self.img, vmin=self.vmin, vmax=self.vmax, aspect='equal', cmap='inferno',alpha=self.mask_alpha) array = self.displayed.get_array().data self.mask = np.zeros(self.img.shape[:2],dtype=bool) @@ -1040,14 +1040,14 @@ class pol_map(object): self.ax = self.fig.add_subplot(111,projection=self.wcs) self.ax_cosmetics() self.cbar_ax = self.fig.add_axes([0.925, 0.12, 0.01, 0.75]) - + #Display selected data (Default to total flux) self.display() #Display polarization vectors in SNR_cut self.pol_vector() #Display integrated values in ROI self.pol_int() - + #Set axes for sliders (SNRp_cut, SNRi_cut) ax_I_cut = self.fig.add_axes([0.125, 0.080, 0.35, 0.01]) ax_P_cut = self.fig.add_axes([0.125, 0.055, 0.35, 0.01]) @@ -1121,7 +1121,7 @@ class pol_map(object): b_crop = Button(ax_crop,"Crop") self.cropped = False b_crop_reset = Button(ax_crop_reset,"Reset") - + def crop(event): if self.cropped: self.cropped = False @@ -1284,7 +1284,7 @@ class pol_map(object): b_snrp.on_clicked(d_snrp) plt.show() - + @property def wcs(self): return deepcopy(WCS(self.Stokes[0].header)) @@ -1312,7 +1312,7 @@ class pol_map(object): @property def data_mask(self): return self.Stokes[np.argmax([self.Stokes[i].header['datatype']=='Data_mask' for i in range(len(self.Stokes))])].data - + def set_data_mask(self, mask): self.Stokes[np.argmax([self.Stokes[i].header['datatype']=='Data_mask' for i in range(len(self.Stokes))])].data = mask.astype(float) @@ -1337,7 +1337,7 @@ class pol_map(object): ax.coords[1].set_axislabel_position('l') ax.coords[1].set_ticklabel_position('l') ax.axis('equal') - + #Display scales and orientation fontprops = fm.FontProperties(size=14) px_size = self.wcs.wcs.cdelt[0]*3600. @@ -1357,7 +1357,7 @@ class pol_map(object): def display(self, fig=None, ax=None): if self.display_selection is None: self.display_selection = "total_flux" - + if self.display_selection.lower() in ['total_flux']: self.data = self.I*self.convert_flux vmin, vmax = 0., np.max(self.data[self.data > 0.]) @@ -1390,12 +1390,12 @@ class pol_map(object): ax = self.ax if hasattr(self, 'im'): self.im.remove() - self.im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno') + self.im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno') self.cbar = plt.colorbar(self.im, cax=self.cbar_ax, label=label) fig.canvas.draw_idle() return self.im else: - im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno') + im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno') ax.set_xlim(0,self.data.shape[1]) ax.set_ylim(0,self.data.shape[0]) plt.colorbar(im, pad=0.025, aspect=80, label=label) @@ -1419,7 +1419,7 @@ class pol_map(object): else: ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=0.5, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.1, color='white') fig.canvas.draw_idle() - + def pol_int(self, fig=None, ax=None): if self.region is None: n_pix = self.I.size @@ -1475,4 +1475,3 @@ class pol_map(object): if not self.region is None: ax.contour(self.region.astype(float),levels=[0.5], colors='white', linewidths=0.8) fig.canvas.draw_idle() - \ No newline at end of file diff --git a/src/lib/reduction.py b/src/lib/reduction.py index bf90217..5d27958 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -482,7 +482,7 @@ def get_error(data_array, headers, error_array=None, data_mask=None, sub_shape=N err_flat = data_array[i]*0.03 error_array[i] = np.sqrt(error_array[i]**2 + error_bkg[i]**2 + err_wav**2 + err_psf**2 + err_flat**2) - + background[i] = sub_image.sum() if (data_array[i] < 0.).any(): print(data_array[i]) @@ -808,7 +808,7 @@ def align_data(data_array, headers, error_array=None, upsample_factor=1., shifts.append(shift) errors.append(error) - + shifts = np.array(shifts) errors = np.array(errors) @@ -818,7 +818,7 @@ def align_data(data_array, headers, error_array=None, upsample_factor=1., for i in range(len(headers_wcs)): headers_wcs[i].wcs.crpix = new_crpix[0] headers[i].update(headers_wcs[i].to_header()) - + data_mask = rescaled_mask.all(axis=0) data_array, error_array, data_mask, headers = crop_array(rescaled_image, headers, rescaled_error, data_mask) @@ -926,7 +926,7 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1., else: raise ValueError("{} is not a valid smoothing option".format(smoothing)) - + return smoothed, error @@ -1039,7 +1039,7 @@ def polarizer_avg(data_array, error_array, data_mask, headers, FWHM=None, err60 = np.sqrt(np.sum(err60_array**2,axis=0)) err120 = np.sqrt(np.sum(err120_array**2,axis=0)) polerr_array = np.array([err0, err60, err120]) - + # Update headers for header in headers: if header['filtnam1']=='POL0': @@ -1186,12 +1186,12 @@ def compute_Stokes(data_array, error_array, data_mask, headers, Q_stokes = np.zeros(pol_array[0].shape) U_stokes = np.zeros(pol_array[0].shape) Stokes_cov = np.zeros((3,3,I_stokes.shape[0],I_stokes.shape[1])) - + for i in range(I_stokes.shape[0]): for j in range(I_stokes.shape[1]): I_stokes[i,j], Q_stokes[i,j], U_stokes[i,j] = np.dot(coeff_stokes, pol_flux[:,i,j]).T Stokes_cov[:,:,i,j] = np.dot(coeff_stokes, np.dot(pol_cov[:,:,i,j], coeff_stokes.T)) - + mask = (Q_stokes**2 + U_stokes**2) > I_stokes**2 if mask.any(): print("WARNING : found {0:d} pixels for which I_pol > I_stokes".format(I_stokes[mask].size)) @@ -1201,17 +1201,17 @@ def compute_Stokes(data_array, error_array, data_mask, headers, dI_dtheta2 = 2.*pol_eff[1]/A*(pol_eff[0]*np.cos(-2.*theta[0]+2.*theta[1])*(pol_flux[2]-I_stokes) - pol_eff[2]*np.cos(-2.*theta[1]+2.*theta[2])*(pol_flux[0]-I_stokes)) dI_dtheta3 = 2.*pol_eff[2]/A*(pol_eff[1]*np.cos(-2.*theta[1]+2.*theta[2])*(pol_flux[0]-I_stokes) - pol_eff[0]*np.cos(-2.*theta[2]+2.*theta[0])*(pol_flux[1]-I_stokes)) dI_dtheta = np.array([dI_dtheta1, dI_dtheta2, dI_dtheta3]) - + dQ_dtheta1 = 2.*pol_eff[0]/A*(np.cos(2.*theta[0])*(pol_flux[1]-pol_flux[2]) - (pol_eff[2]*np.cos(-2.*theta[2]+2.*theta[0]) - pol_eff[1]*np.cos(-2.*theta[0]+2.*theta[1]))*Q_stokes) dQ_dtheta2 = 2.*pol_eff[1]/A*(np.cos(2.*theta[1])*(pol_flux[2]-pol_flux[0]) - (pol_eff[0]*np.cos(-2.*theta[0]+2.*theta[1]) - pol_eff[2]*np.cos(-2.*theta[1]+2.*theta[2]))*Q_stokes) dQ_dtheta3 = 2.*pol_eff[2]/A*(np.cos(2.*theta[2])*(pol_flux[0]-pol_flux[1]) - (pol_eff[1]*np.cos(-2.*theta[1]+2.*theta[2]) - pol_eff[0]*np.cos(-2.*theta[2]+2.*theta[0]))*Q_stokes) dQ_dtheta = np.array([dQ_dtheta1, dQ_dtheta2, dQ_dtheta3]) - + dU_dtheta1 = 2.*pol_eff[0]/A*(np.sin(2.*theta[0])*(pol_flux[1]-pol_flux[2]) - (pol_eff[2]*np.cos(-2.*theta[2]+2.*theta[0]) - pol_eff[1]*np.cos(-2.*theta[0]+2.*theta[1]))*U_stokes) dU_dtheta2 = 2.*pol_eff[1]/A*(np.sin(2.*theta[1])*(pol_flux[2]-pol_flux[0]) - (pol_eff[0]*np.cos(-2.*theta[0]+2.*theta[1]) - pol_eff[2]*np.cos(-2.*theta[1]+2.*theta[2]))*U_stokes) dU_dtheta3 = 2.*pol_eff[2]/A*(np.sin(2.*theta[2])*(pol_flux[0]-pol_flux[1]) - (pol_eff[1]*np.cos(-2.*theta[1]+2.*theta[2]) - pol_eff[0]*np.cos(-2.*theta[2]+2.*theta[0]))*U_stokes) dU_dtheta = np.array([dU_dtheta1, dU_dtheta2, dU_dtheta3]) - + # Compute the uncertainty associated with the polarizers' orientation (see Kishimoto 1999) s_I2_axis = np.sum([dI_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))],axis=0) s_Q2_axis = np.sum([dQ_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))],axis=0) @@ -1221,7 +1221,7 @@ def compute_Stokes(data_array, error_array, data_mask, headers, Stokes_cov[0,0] += s_I2_axis Stokes_cov[1,1] += s_Q2_axis Stokes_cov[2,2] += s_U2_axis - + if not(FWHM is None) and (smoothing.lower() in ['weighted_gaussian_after','weight_gauss_after','gaussian_after','gauss_after']): smoothing = smoothing.lower()[:-6] Stokes_array = np.array([I_stokes, Q_stokes, U_stokes]) @@ -1233,7 +1233,7 @@ 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 - + #Compute integrated values for P, PA before any rotation mask = np.logical_and(data_mask.astype(bool), (I_stokes > 0.)) n_pix = I_stokes[mask].size @@ -1419,11 +1419,11 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, mrot = np.array([[1., 0., 0.], [0., np.cos(2.*alpha), np.sin(2.*alpha)], [0, -np.sin(2.*alpha), np.cos(2.*alpha)]]) - + old_center = np.array(I_stokes.shape)/2 shape = np.fix(np.array(I_stokes.shape)*np.sqrt(2.5)).astype(int) new_center = np.array(shape)/2 - + I_stokes = zeropad(I_stokes, shape) Q_stokes = zeropad(Q_stokes, shape) U_stokes = zeropad(U_stokes, shape) @@ -1433,7 +1433,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, new_Q_stokes = np.zeros(shape) new_U_stokes = np.zeros(shape) new_Stokes_cov = np.zeros((*Stokes_cov.shape[:-2],*shape)) - + for i in range(shape[0]): for j in range(shape[1]): new_I_stokes[i,j], new_Q_stokes[i,j], new_U_stokes[i,j] = np.dot(mrot, np.array([I_stokes[i,j], Q_stokes[i,j], U_stokes[i,j]])).T @@ -1554,7 +1554,7 @@ def rotate_data(data_array, error_array, data_mask, headers, ang): old_center = np.array(data_array[0].shape)/2 shape = np.fix(np.array(data_array[0].shape)*np.sqrt(2.5)).astype(int) new_center = np.array(shape)/2 - + data_array = zeropad(data_array, [data_array.shape[0],*shape]) error_array = zeropad(error_array, [error_array.shape[0],*shape]) data_mask = zeropad(data_mask, shape) diff --git a/src/overplot.py b/src/overplot.py index e695138..11488ce 100755 --- a/src/overplot.py +++ b/src/overplot.py @@ -16,30 +16,30 @@ Stokes_IR = fits.open("../data/IC5063_x3nl030/IR/u2e65g01t_c0f_rot.fits") levelsMorganti = np.array([1.,2.,3.,8.,16.,32.,64.,128.]) -#levels18GHz = np.array([0.6, 1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_18GHz[0].data.max() -levels18GHz = levelsMorganti*0.28*1e-3 -A = overplot_radio(Stokes_UV, Stokes_18GHz) -A.plot(levels=levels18GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/18GHz_overplot_forced.png') +##levels18GHz = np.array([0.6, 1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_18GHz[0].data.max() +#levels18GHz = levelsMorganti*0.28*1e-3 +#A = overplot_radio(Stokes_UV, Stokes_18GHz) +#A.plot(levels=levels18GHz, SNRp_cut=3.0, SNRi_cut=60.0, savename='../plots/IC5063_x3nl030/18GHz_overplot_forced.png') -#levels24GHz = np.array([1.,1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_24GHz[0].data.max() -levels24GHz = levelsMorganti*0.46*1e-3 -B = overplot_radio(Stokes_UV, Stokes_24GHz) -B.plot(levels=levels24GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/24GHz_overplot_forced.png') - -levels103GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_103GHz[0].data[Stokes_103GHz[0].data > 0.])) -C = overplot_radio(Stokes_UV, Stokes_103GHz) -C.plot(levels=levels103GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/103GHz_overplot_forced.png') - -levels229GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_229GHz[0].data[Stokes_229GHz[0].data > 0.])) -D = overplot_radio(Stokes_UV, Stokes_229GHz) -D.plot(levels=levels229GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/229GHz_overplot_forced.png') - -levels357GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_357GHz[0].data[Stokes_357GHz[0].data > 0.])) -E = overplot_radio(Stokes_UV, Stokes_357GHz) -E.plot(levels=levels357GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/357GHz_overplot_forced.png') - -F = overplot_pol(Stokes_UV, Stokes_S2) -F.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/S2_overplot_forced.png', norm=LogNorm(vmin=5e-20,vmax=5e-18)) +##levels24GHz = np.array([1.,1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_24GHz[0].data.max() +#levels24GHz = levelsMorganti*0.46*1e-3 +#B = overplot_radio(Stokes_UV, Stokes_24GHz) +#B.plot(levels=levels24GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/24GHz_overplot_forced.png') +# +#levels103GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_103GHz[0].data[Stokes_103GHz[0].data > 0.])) +#C = overplot_radio(Stokes_UV, Stokes_103GHz) +#C.plot(levels=levels103GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/103GHz_overplot_forced.png') +# +#levels229GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_229GHz[0].data[Stokes_229GHz[0].data > 0.])) +#D = overplot_radio(Stokes_UV, Stokes_229GHz) +#D.plot(levels=levels229GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/229GHz_overplot_forced.png') +# +#levels357GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_357GHz[0].data[Stokes_357GHz[0].data > 0.])) +#E = overplot_radio(Stokes_UV, Stokes_357GHz) +#E.plot(levels=levels357GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/357GHz_overplot_forced.png') +# +#F = overplot_pol(Stokes_UV, Stokes_S2) +#F.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/S2_overplot_forced.png', norm=LogNorm(vmin=5e-20,vmax=5e-18)) G = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r') -G.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/IR_overplot_forced.png', norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r') +G.plot(SNRp_cut=3.0, SNRi_cut=60.0, savename='../plots/IC5063_x3nl030/IR_overplot_forced.png', norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r')