diff --git a/plots/3C405_x136060/3C405_FOC.png b/plots/3C405_x136060/3C405_FOC.png index 52b9d4b..68cf6c0 100644 Binary files a/plots/3C405_x136060/3C405_FOC.png and b/plots/3C405_x136060/3C405_FOC.png differ diff --git a/plots/3C405_x136060/3C405_FOC_SNRi.png b/plots/3C405_x136060/3C405_FOC_SNRi.png new file mode 100644 index 0000000..3b0611e Binary files /dev/null and b/plots/3C405_x136060/3C405_FOC_SNRi.png differ diff --git a/plots/3C405_x136060/3C405_FOC_SNRp.png b/plots/3C405_x136060/3C405_FOC_SNRp.png new file mode 100644 index 0000000..44f06f2 Binary files /dev/null and b/plots/3C405_x136060/3C405_FOC_SNRp.png differ diff --git a/plots/3C405_x136060/3C405_FOC_center_image.png b/plots/3C405_x136060/3C405_FOC_center_image.png index f0c3d60..8ae1106 100644 Binary files a/plots/3C405_x136060/3C405_FOC_center_image.png and b/plots/3C405_x136060/3C405_FOC_center_image.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM070.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM070.png index a1f6987..d8d5570 100644 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM070.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM070.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM070_SNRi.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM070_SNRi.png new file mode 100644 index 0000000..547a629 Binary files /dev/null and b/plots/3C405_x136060/3C405_FOC_combine_FWHM070_SNRi.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM070_SNRp.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM070_SNRp.png new file mode 100644 index 0000000..96862d7 Binary files /dev/null and b/plots/3C405_x136060/3C405_FOC_combine_FWHM070_SNRp.png differ diff --git a/plots/3C405_x136060/3C405_FOC_errors_background_flux.png b/plots/3C405_x136060/3C405_FOC_errors_background_flux.png index 2bde4a0..7d72698 100644 Binary files a/plots/3C405_x136060/3C405_FOC_errors_background_flux.png and b/plots/3C405_x136060/3C405_FOC_errors_background_flux.png differ diff --git a/plots/3C405_x136060/3C405_FOC_errors_background_location.png b/plots/3C405_x136060/3C405_FOC_errors_background_location.png index dea7a7c..a669bc6 100644 Binary files a/plots/3C405_x136060/3C405_FOC_errors_background_location.png and b/plots/3C405_x136060/3C405_FOC_errors_background_location.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_center_image.png b/plots/NGC1068_x274020/NGC1068_FOC_center_image.png index 36421a4..c8cc1c3 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_center_image.png and b/plots/NGC1068_x274020/NGC1068_FOC_center_image.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010.png index 610377b..d432840 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_P.png new file mode 100644 index 0000000..9748cab Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_P_err.png new file mode 100644 index 0000000..0b25070 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_SNRi.png index 9eb1874..f3c7e57 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_SNRi.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_SNRp.png index bb4428d..f1f82ea 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_SNRp.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot.png new file mode 100644 index 0000000..8bbff65 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P.png new file mode 100644 index 0000000..3ef4237 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P_err.png new file mode 100644 index 0000000..9983b33 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRi.png new file mode 100644 index 0000000..32a4802 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRp.png new file mode 100644 index 0000000..d6df816 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2.png index abd815b..dcb19e3 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_P.png new file mode 100644 index 0000000..de72deb Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_P_err.png new file mode 100644 index 0000000..c8136dc Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_SNRi.png index 8bf11f1..4d21573 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_SNRi.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_SNRp.png index c85396e..71375b9 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_SNRp.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot.png new file mode 100644 index 0000000..35cadf1 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot_P.png new file mode 100644 index 0000000..fa7d7c8 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot_P_err.png new file mode 100644 index 0000000..a8f0cc3 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot_SNRi.png new file mode 100644 index 0000000..d04e11c Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot_SNRp.png new file mode 100644 index 0000000..261cf87 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM2_rot_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_errors_background_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_errors_background_flux.png index 8f37c4f..53d7a70 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_errors_background_flux.png and b/plots/NGC1068_x274020/NGC1068_FOC_errors_background_flux.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_errors_background_location.png b/plots/NGC1068_x274020/NGC1068_FOC_errors_background_location.png index 78f052b..f5ead0f 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_errors_background_location.png and b/plots/NGC1068_x274020/NGC1068_FOC_errors_background_location.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010.png new file mode 100644 index 0000000..5ff52b1 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_SNRi.png new file mode 100644 index 0000000..95ca4f1 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_SNRp.png new file mode 100644 index 0000000..78e0cc9 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot.png new file mode 100644 index 0000000..96e9cdb Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot_P.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot_P.png new file mode 100644 index 0000000..263b6ee Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot_P_err.png new file mode 100644 index 0000000..8f4c50b Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot_SNRi.png new file mode 100644 index 0000000..d880d4a Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot_SNRp.png new file mode 100644 index 0000000..9bcb374 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_FWHM010_rot_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010.png new file mode 100644 index 0000000..2c84eb9 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_P.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_P.png new file mode 100644 index 0000000..b38104b Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_P_err.png new file mode 100644 index 0000000..11d6618 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_SNRi.png new file mode 100644 index 0000000..76c3f5b Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_SNRp.png new file mode 100644 index 0000000..a32010a Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot.png new file mode 100644 index 0000000..abfdc8d Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot_P.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot_P.png new file mode 100644 index 0000000..8c5f064 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot_P_err.png new file mode 100644 index 0000000..744d318 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot_SNRi.png new file mode 100644 index 0000000..05a645f Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot_SNRp.png new file mode 100644 index 0000000..d7802f5 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_gaussian_after_FWHM010_rot_SNRp.png differ diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 4ede500..6000729 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -77,22 +77,21 @@ def main(): # Data binning rebin = True if rebin: - pxsize = 0.1 + pxsize = 0.10 px_scale = 'arcsec' #pixel or arcsec rebin_operation = 'sum' #sum or average # Alignement align_center = 'image' #If None will align image to image center display_data = False # Smoothing - smoothing_function = 'combine' #gaussian or combine - smoothing_FWHM = 0.1 #If None, no smoothing is done + smoothing_function = 'gaussian_after' #gaussian_after, gaussian or combine + smoothing_FWHM = 0.10 #If None, no smoothing is done smoothing_scale = 'arcsec' #pixel or arcsec # Rotation rotate = False #rotation to North convention can give erroneous results - rotate_library = 'scipy' #scipy or pillow # Polarization map output figname = 'NGC1068_FOC' #target/intrument name - figtype = '_combine_FWHM010' #additionnal informations + figtype = '_gaussian_after_FWHM010' #additionnal informations SNRp_cut = 3 #P measurments with SNR>3 SNRi_cut = 30 #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 @@ -130,10 +129,7 @@ def main(): # Rotate images to have North up if rotate: ref_header = copy.deepcopy(headers[0]) - if rotate_library.lower() in ['scipy','scipy.ndimage']: - I_stokes, Q_stokes, U_stokes, Stokes_cov, headers = proj_red.rotate_sc_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, -ref_header['orientat']) - if rotate_library.lower() in ['pillow','pil']: - I_stokes, Q_stokes, U_stokes, Stokes_cov, headers = proj_red.rotate_PIL_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, -ref_header['orientat']) + I_stokes, Q_stokes, U_stokes, Stokes_cov, headers = proj_red.rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, -ref_header['orientat']) # 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) @@ -144,12 +140,12 @@ def main(): ## Step 5: # Plot polarization map (Background is either total Flux, Polarization degree or Polarization degree error). proj_plots.polarization_map(copy.deepcopy(Stokes_test), SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype, plots_folder=plots_folder, display=None) - #proj_plots.polarization_map(copy.deepcopy(Stokes_test), SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P", plots_folder=plots_folder, display='Pol_deg') - #proj_plots.polarization_map(copy.deepcopy(Stokes_test), SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P_err", plots_folder=plots_folder, display='Pol_deg_err') + proj_plots.polarization_map(copy.deepcopy(Stokes_test), SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P", plots_folder=plots_folder, display='Pol_deg') + proj_plots.polarization_map(copy.deepcopy(Stokes_test), SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P_err", plots_folder=plots_folder, display='Pol_deg_err') proj_plots.polarization_map(copy.deepcopy(Stokes_test), SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_SNRi", plots_folder=plots_folder, display='SNRi') proj_plots.polarization_map(copy.deepcopy(Stokes_test), SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_SNRp", plots_folder=plots_folder, display='SNRp') return 0 -#if __name__ == "__main__": -# sys.exit(main()) +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/FOC_reduction2_ELR.py b/src/FOC_reduction2_ELR.py index 20a8270..2d80c97 100755 --- a/src/FOC_reduction2_ELR.py +++ b/src/FOC_reduction2_ELR.py @@ -16,11 +16,11 @@ import lib.plots_ELR as proj_plots #Functions for plotting data 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'] - 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'] +# globals()['plots_folder'] = "../plots/NGC1068_x274020/" # globals()['data_folder'] = "../data/NGC1068_x14w010/" # infiles = ['x14w0101t_c0f.fits','x14w0102t_c0f.fits','x14w0103t_c0f.fits', @@ -29,10 +29,10 @@ def main(): # 'x14w0104t_c1f.fits','x14w0105p_c1f.fits','x14w0106t_c1f.fits'] # globals()['plots_folder'] = "../plots/NGC1068_x14w010/" -# globals()['data_folder'] = "../data/3C405_x136060/" -# infiles = ['x1360601t_c0f.fits','x1360602t_c0f.fits','x1360603t_c0f.fits'] + globals()['data_folder'] = "../data/3C405_x136060/" + infiles = ['x1360601t_c0f.fits','x1360602t_c0f.fits','x1360603t_c0f.fits'] # infiles = ['x1360601t_c1f.fits','x1360602t_c1f.fits','x1360603t_c1f.fits'] -# globals()['plots_folder'] = "../plots/3C405_x136060/" + globals()['plots_folder'] = "../plots/3C405_x136060/" # globals()['data_folder'] = "../data/CygnusA_x43w0/" # infiles = ['x43w0101r_c0f.fits', 'x43w0104r_c0f.fits', 'x43w0107r_c0f.fits', @@ -72,17 +72,17 @@ def main(): psf_shape=(9,9) iterations = 10 # Error estimation - error_sub_shape = (75,75) - display_error = False + error_sub_shape = (150,150) + display_error = True # Data binning rebin = True if rebin: - pxsize = 8 + pxsize = 15 px_scale = 'pixel' #pixel or arcsec rebin_operation = 'average' #sum or average # Alignement align_center = 'image' #If None will align image to image center - display_data = False + display_data = True # Smoothing smoothing_function = 'gaussian' #gaussian or combine smoothing_FWHM = 1 #If None, no smoothing is done @@ -91,8 +91,8 @@ def main(): rotate = False #rotation to North convention can give erroneous results rotate_library = 'scipy' #scipy or pillow # Polarization map output - figname = 'NGC1068_FOC' #target/intrument name - figtype = '_ELR_same_param_same_op_sP100' #additionnal informations + figname = '3C405_FOC' #target/intrument name + figtype = '' #additionnal informations SNRp_cut = 3 #P measurments with SNR>3 SNRi_cut = 4 #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 diff --git a/src/lib/fits.py b/src/lib/fits.py index 32dedf0..f4d6fed 100755 --- a/src/lib/fits.py +++ b/src/lib/fits.py @@ -43,8 +43,8 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): if compute_flux: for i in range(len(infiles)): - # Compute the flux in ergs/cm²/s.angstrom - data_array[i] *= headers[i]['PHOTFLAM']/headers[i]['EXPTIME'] + # Compute the flux in counts/sec + data_array[i] /= headers[i]['EXPTIME'] return data_array, headers diff --git a/src/lib/plots.py b/src/lib/plots.py index 40490a7..43b315a 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -134,12 +134,12 @@ def polarization_map(Stokes, SNRp_cut=3., SNRi_cut=30., step_vec=1, pang_err = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_ang_err' for i in range(len(Stokes))])] pivot_wav = Stokes[0].header['photplam'] + convert_flux = Stokes[0].header['photflam'] wcs = WCS(Stokes[0]).deepcopy() #Compute SNR and apply cuts pol.data[pol.data == 0.] = np.nan SNRp = pol.data/pol_err.data - #SNRp = pol.data/pol_err_Poisson.data pol.data[SNRp < SNRp_cut] = np.nan SNRi = stkI.data/np.sqrt(stk_cov.data[0,0]) pol.data[SNRi < SNRi_cut] = np.nan @@ -162,8 +162,8 @@ def polarization_map(Stokes, SNRp_cut=3., SNRi_cut=30., step_vec=1, if display is None: # If no display selected, show intensity map - vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]) - im = ax.imshow(stkI.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.) + vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux) + im = ax.imshow(stkI.data*convert_flux,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}$ [$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) @@ -193,11 +193,11 @@ def polarization_map(Stokes, SNRp_cut=3., SNRi_cut=30., step_vec=1, cont = ax.contour(SNRp, extent=[-SNRp.shape[1]/2.,SNRp.shape[1]/2.,-SNRp.shape[0]/2.,SNRp.shape[0]/2.], levels=levelsP, colors='grey', linewidths=0.5) else: # Defaults to intensity map - vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]) - im = ax.imshow(stkI.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.) + vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux) + im = ax.imshow(stkI.data*convert_flux,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}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") levelsI = np.linspace(SNRi_cut, SNRi.max(), 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) - cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") px_size = wcs.wcs.get_cdelt()[0] px_sc = AnchoredSizeBar(ax.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w') @@ -226,7 +226,7 @@ def polarization_map(Stokes, SNRp_cut=3., SNRi_cut=30., step_vec=1, PA_int = (90./np.pi)*np.arctan2(U_int,Q_int)+90. 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) - ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1:.1e} $\pm$ {2:.1e} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(pivot_wav,I_int,I_int_err)+"\n"+r"$P^{{int}}$ = {0:.2f} $\pm$ {1:.2f} %".format(P_int,P_int_err)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.2f} $\pm$ {1:.2f} °".format(PA_int,PA_int_err), color='white', fontsize=11, xy=(0.01, 0.94), xycoords='axes fraction') + ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1:.1e} $\pm$ {2:.1e} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(pivot_wav,I_int*convert_flux,I_int_err*convert_flux)+"\n"+r"$P^{{int}}$ = {0:.2f} $\pm$ {1:.2f} %".format(P_int,P_int_err)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.2f} $\pm$ {1:.2f} °".format(PA_int,PA_int_err), color='white', fontsize=11, xy=(0.01, 0.94), 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/plots_ELR.py b/src/lib/plots_ELR.py index cff5904..2f69476 100755 --- a/src/lib/plots_ELR.py +++ b/src/lib/plots_ELR.py @@ -196,7 +196,7 @@ def polarization_map(Stokes, SNRp_cut=3., SNRi_cut=30., step_vec=1, else: # Defaults to intensity map vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux) - im = ax.imshow(stkI.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.) + im = ax.imshow(stkI.data*convert_flux,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}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") levelsI = np.linspace(SNRi_cut, SNRi.max(), 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) diff --git a/src/lib/reduction.py b/src/lib/reduction.py index a234696..e7f1aa4 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -33,19 +33,12 @@ prototypes : Compute polarization degree (in %) and angle (in degree) and their respective errors - - rotate_sc_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang) -> I_stokes, Q_stokes, U_stokes, Stokes_cov, headers + - rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang) -> I_stokes, Q_stokes, U_stokes, Stokes_cov, headers Rotate I, Q, U given an angle in degrees using scipy functions. - - rotate_PIL_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang) -> I_stokes, Q_stokes, U_stokes, Stokes_cov, headers - Rotate I, Q, U given an angle in degrees using Pillow function. - - - rotate2_sc_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang) -> I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, headers + - rotate2_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang) -> I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, headers Rotate I, Q, U, P, PA and associated errors given an angle in degrees using scipy functions. - - - rotate2_PIL_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang) -> I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, headers - Rotate I, Q, U, P, PA and associated errors given an angle in degrees - using Pillow function. """ import copy @@ -55,7 +48,6 @@ import matplotlib.dates as mdates from datetime import datetime from scipy.ndimage import rotate as sc_rotate from scipy.ndimage import shift as sc_shift -from PIL import Image from astropy.wcs import WCS from lib.deconvolve import deconvolve_im, gaussian_psf from lib.convex_hull import image_hull @@ -63,36 +55,6 @@ from lib.plots import plot_obs from lib.cross_correlation import phase_cross_correlation -def PIL_rotate(ndarray, angle, reshape=False, cval=0.): - """ - Define Function with similar parameters as scipy.ndimage.rotate making use - of the Pillow library. - ---------- - Inputs: - ndarray : numpy.ndarray - 2D float array to be rotated using Pillow library functions. - angle : float - Counter-clockwise ngle in degrees that the input array should be - rotated to. - reshape : boolean, optional - If True, output image will be of shape such that the full rotated - input array is displayed. If False, output array will be of same - shape than input array, cutting rotated edges. - Defaults to False. - cval : float, optional - Values with which will be filled pixels entering the output array by - rotation of the input array. - Dafaults to 0. - ---------- - Returns: - rotated_array : numpy.ndarray - Rotated array using Pillow functions. - """ - image = Image.fromarray(ndarray) - rotated = image.rotate(angle, expand=reshape, fillcolor=cval) - return np.asarray(rotated) - - def get_row_compressor(old_dimension, new_dimension, operation='sum'): """ Return the matrix that allows to compress an array from an old dimension of @@ -660,8 +622,8 @@ def align_data(data_array, error_array=None, upsample_factor=1., ref_data=None, full_array = np.concatenate((data_array,[ref_data]),axis=0) err_array = np.concatenate((error_array,[np.zeros(ref_data.shape)]),axis=0) - #full_array, err_array = crop_array(full_array, err_array, step=5, - # inside=False) + full_array, err_array = crop_array(full_array, err_array, step=5, + inside=False) data_array, ref_data = full_array[:-1], full_array[-1] error_array = err_array[:-1] @@ -772,7 +734,7 @@ def smooth_data(data_array, error_array, headers, FWHM=1., scale='pixel', px_dim[0] = 50. w.wcs.cdelt = 206.3*px_dim/(f_ratio*HST_aper) header.update(w.to_header()) - pxsize[i] = np.round(w.wcs.cdelt,5) + pxsize[i] = np.round(w.wcs.cdelt,4) if (pxsize != pxsize[0]).any(): raise ValueError("Not all images in array have same pixel size") FWHM /= pxsize[0].min() @@ -915,8 +877,17 @@ def polarizer_avg(data_array, error_array, headers, FWHM=None, scale='pixel', err120 = np.mean(err120_array,axis=0)/np.sqrt(err120_array.shape[0]) polerr_array = np.array([err0, err60, err120]) - headers_array = headers0 + headers60 + headers120 - if not(FWHM is None): + # Update headers + for header in headers: + if header['filtnam1']=='POL0': + list_head = headers0 + elif header['filtnam1']=='POL60': + list_head = headers60 + else: + list_head = headers120 + header['exptime'] = np.mean([head['exptime'] for head in list_head])/len(list_head) + headers_array = [headers0[0], headers60[0], headers120[0]] + if not(FWHM is None) and (smoothing.lower() in ['gaussian','gauss']): # Smooth by convoluting with a gaussian each polX image. pol_array, polerr_array = smooth_data(pol_array, polerr_array, headers_array, FWHM=FWHM, scale=scale) @@ -943,7 +914,7 @@ def polarizer_avg(data_array, error_array, headers, FWHM=None, scale='pixel', def compute_Stokes(data_array, error_array, headers, FWHM=None, - scale='pixel', smoothing='gaussian'): + scale='pixel', smoothing='gaussian_after'): """ Compute the Stokes parameters I, Q and U for a given data_set ---------- @@ -970,7 +941,9 @@ def compute_Stokes(data_array, error_array, headers, FWHM=None, weighted pixels (inverse square of errors). -'gaussian' convolve any input image with a gaussian of standard deviation stdev = FWHM/(2*sqrt(2*log(2))). - Defaults to 'gaussian'. Won't be used if FWHM is None. + -'gaussian_after' convolve output Stokes I/Q/U with a gaussian of + standard deviation stdev = FWHM/(2*sqrt(2*log(2))). + Defaults to 'gaussian_after'. Won't be used if FWHM is None. ---------- Returns: I_stokes : numpy.ndarray @@ -1021,6 +994,17 @@ def compute_Stokes(data_array, error_array, headers, FWHM=None, Q_stokes[np.isnan(Q_stokes)]=0. U_stokes[np.isnan(U_stokes)]=0. + if not(FWHM is None) and (smoothing.lower() in ['gaussian_after','gauss_after']): + Stokes_array = np.array([I_stokes, Q_stokes, U_stokes]) + Stokes_error = np.array([np.sqrt(Stokes_cov[i,i]) for i in range(3)]) + Stokes_headers = headers[0:3] + + Stokes_array, Stokes_error = smooth_data(Stokes_array, Stokes_error, + headers=Stokes_headers, FWHM=FWHM, scale=scale) + + 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 @@ -1092,7 +1076,7 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers): return P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P -def rotate_sc_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang): +def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang): """ 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 @@ -1184,99 +1168,7 @@ def rotate_sc_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang): return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_headers -def rotate_PIL_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang): - """ - 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 - orientation keyword. - ---------- - Inputs: - I_stokes : numpy.ndarray - Image (2D floats) containing the Stokes parameters accounting for - total intensity - Q_stokes : numpy.ndarray - Image (2D floats) containing the Stokes parameters accounting for - vertical/horizontal linear polarization intensity - U_stokes : numpy.ndarray - Image (2D floats) containing the Stokes parameters accounting for - +45/-45deg linear polarization intensity - Stokes_cov : numpy.ndarray - Covariance matrix of the Stokes parameters I, Q, U. - headers : header list - List of headers corresponding to the reduced images. - ang : float - Rotation angle (in degrees) that should be applied to the Stokes - parameters - ---------- - Returns: - new_I_stokes : numpy.ndarray - Rotated mage (2D floats) containing the rotated Stokes parameters - accounting for total intensity - new_Q_stokes : numpy.ndarray - Rotated mage (2D floats) containing the rotated Stokes parameters - accounting for vertical/horizontal linear polarization intensity - new_U_stokes : numpy.ndarray - Rotated image (2D floats) containing the rotated Stokes parameters - accounting for +45/-45deg linear polarization intensity. - new_Stokes_cov : numpy.ndarray - Updated covariance matrix of the Stokes parameters I, Q, U. - new_headers : header list - Updated list of headers corresponding to the reduced images accounting - for the new orientation angle. - """ - #Rotate I_stokes, Q_stokes, U_stokes using rotation matrix - alpha = ang*np.pi/180. - new_I_stokes = 1.*I_stokes - 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 - - #Compute new covariance matrix on rotated parameters - new_Stokes_cov = copy.deepcopy(Stokes_cov) - new_Stokes_cov[1,1] = np.cos(2.*alpha)**2*Stokes_cov[1,1] + np.sin(2.*alpha)**2*Stokes_cov[2,2] + 2.*np.cos(2.*alpha)*np.sin(2.*alpha)*Stokes_cov[1,2] - new_Stokes_cov[2,2] = np.sin(2.*alpha)**2*Stokes_cov[1,1] + np.cos(2.*alpha)**2*Stokes_cov[2,2] - 2.*np.cos(2.*alpha)*np.sin(2.*alpha)*Stokes_cov[1,2] - new_Stokes_cov[0,1] = new_Stokes_cov[1,0] = np.cos(2.*alpha)*Stokes_cov[0,1] + np.sin(2.*alpha)*Stokes_cov[0,2] - new_Stokes_cov[0,2] = new_Stokes_cov[2,0] = -np.sin(2.*alpha)*Stokes_cov[0,1] + np.cos(2.*alpha)*Stokes_cov[0,2] - new_Stokes_cov[1,2] = new_Stokes_cov[2,1] = np.cos(2.*alpha)*np.sin(2.*alpha)*(Stokes_cov[2,2] - Stokes_cov[1,1]) + (np.cos(2.*alpha)**2 - np.sin(2.*alpha)**2)*Stokes_cov[1,2] - - #Rotate original images using scipy.ndimage.rotate - new_I_stokes = PIL_rotate(new_I_stokes, ang, reshape=False, - cval=np.sqrt(new_Stokes_cov[0,0][0,0])) - new_Q_stokes = PIL_rotate(new_Q_stokes, ang, reshape=False, - cval=np.sqrt(new_Stokes_cov[1,1][0,0])) - new_U_stokes = PIL_rotate(new_U_stokes, ang, reshape=False, - cval=np.sqrt(new_Stokes_cov[2,2][0,0])) - for i in range(3): - for j in range(3): - new_Stokes_cov[i,j] = PIL_rotate(new_Stokes_cov[i,j], ang, - reshape=False, cval=new_Stokes_cov[i,j].mean()) - - #Update headers to new angle - new_headers = [] - mrot = np.array([[np.cos(-alpha), -np.sin(-alpha)], - [np.sin(-alpha), np.cos(-alpha)]]) - for header in headers: - new_header = copy.deepcopy(header) - new_header['orientat'] = header['orientat'] + ang - - new_wcs = WCS(header).deepcopy() - if new_wcs.wcs.has_cd(): # CD matrix - del w.wcs.cd - keys = ['CD1_1','CD1_2','CD2_1','CD2_2'] - for key in keys: - new_header.remove(key, ignore_missing=True) - w.wcs.cdelt = 3600.*np.sqrt(np.sum(w.wcs.get_pc()**2,axis=1)) - elif new_wcs.wcs.has_pc(): # PC matrix + CDELT - newpc = np.dot(mrot, new_wcs.wcs.get_pc()) - new_wcs.wcs.pc = newpc - new_wcs.wcs.set() - new_header.update(new_wcs.to_header()) - - new_headers.append(new_header) - - return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_headers - - -def rotate2_sc_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang): +def rotate2_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang): """ 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 @@ -1391,120 +1283,3 @@ def rotate2_sc_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang): new_headers.append(new_header) return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, new_headers - - -def rotate2_PIL_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, ang): - """ - 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 - orientation keyword. - ---------- - Inputs: - I_stokes : numpy.ndarray - Image (2D floats) containing the Stokes parameters accounting for - total intensity - Q_stokes : numpy.ndarray - Image (2D floats) containing the Stokes parameters accounting for - vertical/horizontal linear polarization intensity - U_stokes : numpy.ndarray - Image (2D floats) containing the Stokes parameters accounting for - +45/-45deg linear polarization intensity - Stokes_cov : numpy.ndarray - Covariance matrix of the Stokes parameters I, Q, U. - headers : header list - List of headers corresponding to the reduced images. - ang : float - Rotation angle (in degrees) that should be applied to the Stokes - parameters - ---------- - Returns: - new_I_stokes : numpy.ndarray - Rotated mage (2D floats) containing the rotated Stokes parameters - accounting for total intensity - new_Q_stokes : numpy.ndarray - Rotated mage (2D floats) containing the rotated Stokes parameters - accounting for vertical/horizontal linear polarization intensity - new_U_stokes : numpy.ndarray - Rotated image (2D floats) containing the rotated Stokes parameters - accounting for +45/-45deg linear polarization intensity. - new_Stokes_cov : numpy.ndarray - Updated covariance matrix of the Stokes parameters I, Q, U. - P : numpy.ndarray - Image (2D floats) containing the polarization degree (in %). - s_P : numpy.ndarray - Image (2D floats) containing the error on the polarization degree. - PA : numpy.ndarray - Image (2D floats) containing the polarization angle. - s_PA : numpy.ndarray - Image (2D floats) containing the error on the polarization angle. - debiased_P : numpy.ndarray - Image (2D floats) containing the debiased polarization degree (in %). - s_P_P : numpy.ndarray - Image (2D floats) containing the Poisson noise error on the - polarization degree. - s_PA_P : numpy.ndarray - Image (2D floats) containing the Poisson noise error on the - polarization angle. - """ - # Rotate I_stokes, Q_stokes, U_stokes using rotation matrix - alpha = ang*np.pi/180. - new_I_stokes = 1.*I_stokes - 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 - - # Compute new covariance matrix on rotated parameters - new_Stokes_cov = copy.deepcopy(Stokes_cov) - new_Stokes_cov[1,1] = np.cos(2.*alpha)**2*Stokes_cov[1,1] + np.sin(2.*alpha)**2*Stokes_cov[2,2] + 2.*np.cos(2.*alpha)*np.sin(2.*alpha)*Stokes_cov[1,2] - new_Stokes_cov[2,2] = np.sin(2.*alpha)**2*Stokes_cov[1,1] + np.cos(2.*alpha)**2*Stokes_cov[2,2] - 2.*np.cos(2.*alpha)*np.sin(2.*alpha)*Stokes_cov[1,2] - new_Stokes_cov[0,1] = new_Stokes_cov[1,0] = np.cos(2.*alpha)*Stokes_cov[0,1] + np.sin(2.*alpha)*Stokes_cov[0,2] - new_Stokes_cov[0,2] = new_Stokes_cov[2,0] = -np.sin(2.*alpha)*Stokes_cov[0,1] + np.cos(2.*alpha)*Stokes_cov[0,2] - new_Stokes_cov[1,2] = new_Stokes_cov[2,1] = np.cos(2.*alpha)*np.sin(2.*alpha)*(Stokes_cov[2,2] - Stokes_cov[1,1]) + (np.cos(2.*alpha)**2 - np.sin(2.*alpha)**2)*Stokes_cov[1,2] - - # Compute new polarization parameters - P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = compute_pol(new_I_stokes, - new_Q_stokes, new_U_stokes, new_Stokes_cov, headers) - - # Rotate original images using scipy.ndimage.rotate - new_I_stokes = PIL_rotate(new_I_stokes, ang, reshape=False, - cval=np.sqrt(new_Stokes_cov[0,0][0,0])) - new_Q_stokes = PIL_rotate(new_Q_stokes, ang, reshape=False, - cval=np.sqrt(new_Stokes_cov[1,1][0,0])) - new_U_stokes = PIL_rotate(new_U_stokes, ang, reshape=False, - cval=np.sqrt(new_Stokes_cov[2,2][0,0])) - P = PIL_rotate(P, ang, reshape=False, cval=P.mean())+0. - debiased_P = PIL_rotate(debiased_P, ang, reshape=False, - cval=debiased_P.mean()) - s_P = PIL_rotate(s_P, ang, reshape=False, cval=s_P.mean()) - s_P_P = PIL_rotate(s_P_P, ang, reshape=False, cval=s_P_P.mean()) - PA = PIL_rotate(PA, ang, reshape=False, cval=PA.mean()) - s_PA = PIL_rotate(s_PA, ang, reshape=False, cval=s_PA.mean()) - s_PA_P = PIL_rotate(s_PA_P, ang, reshape=False, cval=s_PA_P.mean()) - for i in range(3): - for j in range(3): - new_Stokes_cov[i,j] = PIL_rotate(new_Stokes_cov[i,j], ang, - reshape=False, cval=new_Stokes_cov[i,j].mean()) - - #Update headers to new angle - new_headers = [] - mrot = np.array([[np.cos(-alpha), -np.sin(-alpha)], - [np.sin(-alpha), np.cos(-alpha)]]) - for header in headers: - new_header = copy.deepcopy(header) - new_header['orientat'] = header['orientat'] + ang - - new_wcs = WCS(header).deepcopy() - if new_wcs.wcs.has_cd(): # CD matrix - del w.wcs.cd - keys = ['CD1_1','CD1_2','CD2_1','CD2_2'] - for key in keys: - new_header.remove(key, ignore_missing=True) - w.wcs.cdelt = 3600.*np.sqrt(np.sum(w.wcs.get_pc()**2,axis=1)) - elif new_wcs.wcs.has_pc(): # PC matrix + CDELT - newpc = np.dot(mrot, new_wcs.wcs.get_pc()) - new_wcs.wcs.pc = newpc - new_wcs.wcs.set() - new_header.update(new_wcs.to_header()) - - new_headers.append(new_header) - - return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P, new_headers