diff --git a/plots/IC5063_x3nl030/18GHz_overplot.png b/plots/IC5063_x3nl030/18GHz_overplot.png index 7b95a91..ef961e4 100644 Binary files a/plots/IC5063_x3nl030/18GHz_overplot.png and b/plots/IC5063_x3nl030/18GHz_overplot.png differ diff --git a/plots/IC5063_x3nl030/18GHz_overplot_forced.png b/plots/IC5063_x3nl030/18GHz_overplot_forced.png index 3e92a92..5ebd332 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/18GHz_overplot_forced_maxUV.png b/plots/IC5063_x3nl030/18GHz_overplot_forced_maxUV.png index 1dabbef..40ff2dc 100644 Binary files a/plots/IC5063_x3nl030/18GHz_overplot_forced_maxUV.png and b/plots/IC5063_x3nl030/18GHz_overplot_forced_maxUV.png differ diff --git a/plots/IC5063_x3nl030/24GHz_overplot.png b/plots/IC5063_x3nl030/24GHz_overplot.png index c8aeeb2..9a02140 100644 Binary files a/plots/IC5063_x3nl030/24GHz_overplot.png and b/plots/IC5063_x3nl030/24GHz_overplot.png differ diff --git a/plots/IC5063_x3nl030/24GHz_overplot_forced.png b/plots/IC5063_x3nl030/24GHz_overplot_forced.png index c82185d..d4da6e8 100644 Binary files a/plots/IC5063_x3nl030/24GHz_overplot_forced.png and b/plots/IC5063_x3nl030/24GHz_overplot_forced.png differ diff --git a/plots/IC5063_x3nl030/24GHz_overplot_forced_maxUV.png b/plots/IC5063_x3nl030/24GHz_overplot_forced_maxUV.png index a759cb3..933d632 100644 Binary files a/plots/IC5063_x3nl030/24GHz_overplot_forced_maxUV.png and b/plots/IC5063_x3nl030/24GHz_overplot_forced_maxUV.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_1px.png b/plots/IC5063_x3nl030/IC5063_FOC_1px.png new file mode 100644 index 0000000..2120464 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_1px.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_1px_IQU.png b/plots/IC5063_x3nl030/IC5063_FOC_1px_IQU.png new file mode 100644 index 0000000..c8dc758 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_1px_IQU.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_POL0_crop_region.png b/plots/IC5063_x3nl030/IC5063_FOC_POL0_crop_region.png index 10176ca..524ed95 100755 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_POL0_crop_region.png and b/plots/IC5063_x3nl030/IC5063_FOC_POL0_crop_region.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_center_image.png b/plots/IC5063_x3nl030/IC5063_FOC_center_image.png index f37dc85..b0b2c45 100755 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_center_image.png and b/plots/IC5063_x3nl030/IC5063_FOC_center_image.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae.png new file mode 100644 index 0000000..a929945 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_IQU.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_IQU.png new file mode 100644 index 0000000..58f25da Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_IQU.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_I_err.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_I_err.png new file mode 100644 index 0000000..7b95d3f Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_I_err.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_P.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_P.png new file mode 100644 index 0000000..bb56063 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_P.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_P_err.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_P_err.png new file mode 100644 index 0000000..60db942 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_P_err.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_P_flux.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_P_flux.png new file mode 100644 index 0000000..4a433ac Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_P_flux.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_SNRi.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_SNRi.png new file mode 100644 index 0000000..b37dc11 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_SNRi.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_SNRp.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_SNRp.png new file mode 100644 index 0000000..0857b46 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_wae_SNRp.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_crop_region.png b/plots/IC5063_x3nl030/IC5063_FOC_crop_region.png index d6b9115..3e443a9 100755 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_crop_region.png and b/plots/IC5063_x3nl030/IC5063_FOC_crop_region.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_errors_background_flux.png b/plots/IC5063_x3nl030/IC5063_FOC_errors_background_flux.png index 2d9e4e1..6ec65c3 100755 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_errors_background_flux.png and b/plots/IC5063_x3nl030/IC5063_FOC_errors_background_flux.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_errors_background_location.png b/plots/IC5063_x3nl030/IC5063_FOC_errors_background_location.png index 8bf2812..519dc76 100755 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_errors_background_location.png and b/plots/IC5063_x3nl030/IC5063_FOC_errors_background_location.png differ diff --git a/plots/NGC1068_x274020/Figure_1.png b/plots/NGC1068_x274020/Figure_1.png deleted file mode 100644 index 35523d4..0000000 Binary files a/plots/NGC1068_x274020/Figure_1.png and /dev/null differ diff --git a/plots/NGC1068_x274020/Figure_2.png b/plots/NGC1068_x274020/Figure_2.png deleted file mode 100644 index 3c56a0f..0000000 Binary files a/plots/NGC1068_x274020/Figure_2.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC.png b/plots/NGC1068_x274020/NGC1068_FOC.png index 6eaa738..794bc22 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC.png and b/plots/NGC1068_x274020/NGC1068_FOC.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_IQU.png b/plots/NGC1068_x274020/NGC1068_FOC_IQU.png index a4d6c01..7963080 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_IQU.png and b/plots/NGC1068_x274020/NGC1068_FOC_IQU.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_I_err.png b/plots/NGC1068_x274020/NGC1068_FOC_I_err.png new file mode 100644 index 0000000..3c482e9 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_I_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_P.png b/plots/NGC1068_x274020/NGC1068_FOC_P.png index bc17214..8c81abd 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_P.png and b/plots/NGC1068_x274020/NGC1068_FOC_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_POL0_crop_region.png b/plots/NGC1068_x274020/NGC1068_FOC_POL0_crop_region.png index 7123000..cc248e5 100755 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_POL0_crop_region.png and b/plots/NGC1068_x274020/NGC1068_FOC_POL0_crop_region.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_P_err.png index 145a4cc..24e236e 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_P_err.png and b/plots/NGC1068_x274020/NGC1068_FOC_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_P_flux.png index c3ec372..b9a0a51 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_P_flux.png and b/plots/NGC1068_x274020/NGC1068_FOC_P_flux.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_SNRi.png index c39e530..f65c1be 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_SNRi.png and b/plots/NGC1068_x274020/NGC1068_FOC_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_SNRp.png index 53225b8..dfc2622 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_SNRp.png and b/plots/NGC1068_x274020/NGC1068_FOC_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_center_image.png b/plots/NGC1068_x274020/NGC1068_FOC_center_image.png index d83bb2c..432e75c 100755 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_FWHM005.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005.png new file mode 100644 index 0000000..220d740 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_IQU.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_IQU.png new file mode 100644 index 0000000..64cc509 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_IQU.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_I_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_I_err.png new file mode 100644 index 0000000..7b7b613 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_I_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_P.png new file mode 100644 index 0000000..ab2a93d Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_P_err.png new file mode 100644 index 0000000..321a488 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_P_flux.png new file mode 100644 index 0000000..4ffbec1 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_P_flux.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_SNRi.png new file mode 100644 index 0000000..618f6a3 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_SNRp.png new file mode 100644 index 0000000..19e8afb Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM005_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020.png index 61ce1eb..d627b97 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_I_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_I_err.png index 6a0da42..22bc74f 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_I_err.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_I_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P.png index b8f657b..1ca5ba7 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_err.png index 0d7f57c..18d226f 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_err.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_flux.png index 58d0590..7148bdf 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_flux.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_flux.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRi.png index 6395b3f..adf3902 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRi.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRp.png index c5aa6c8..ce31bfe 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRp.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test.png deleted file mode 100644 index ed80390..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_IQU.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_IQU.png deleted file mode 100644 index 30d049b..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_IQU.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_I_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_I_err.png deleted file mode 100644 index 9c12ce4..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_I_err.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_P.png deleted file mode 100644 index 21c2353..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_P.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_P_err.png deleted file mode 100644 index ce4e805..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_P_err.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_P_flux.png deleted file mode 100644 index 72425fa..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_P_flux.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_SNRi.png deleted file mode 100644 index 70ae5c0..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_SNRi.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_SNRp.png deleted file mode 100644 index 8a4045a..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_SNRp.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae.png deleted file mode 100644 index 2811e9b..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_IQU.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_IQU.png deleted file mode 100644 index a815c09..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_IQU.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_I_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_I_err.png deleted file mode 100644 index aacec9f..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_I_err.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_P.png deleted file mode 100644 index 5b12718..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_P.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_P_err.png deleted file mode 100644 index 7721803..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_P_err.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_P_flux.png deleted file mode 100644 index 9b7af30..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_P_flux.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_SNRi.png deleted file mode 100644 index 3a66585..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_SNRi.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_SNRp.png deleted file mode 100644 index ad50aa7..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_test_wae_SNRp.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae.png new file mode 100644 index 0000000..4a49a72 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP.png new file mode 100644 index 0000000..6bdcf9e Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_IQU.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_IQU.png new file mode 100644 index 0000000..60d2ba5 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_IQU.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_I_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_I_err.png new file mode 100644 index 0000000..58abd9c Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_I_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_P.png new file mode 100644 index 0000000..bd4c277 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_P_err.png new file mode 100644 index 0000000..9ef22c7 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_P_flux.png new file mode 100644 index 0000000..d1d9073 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_P_flux.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_SNRi.png new file mode 100644 index 0000000..75e6c97 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_SNRp.png new file mode 100644 index 0000000..9543c49 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_waeP_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_IQU.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_IQU.png new file mode 100644 index 0000000..60d2ba5 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_IQU.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_I_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_I_err.png new file mode 100644 index 0000000..47f286f Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_I_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_P.png new file mode 100644 index 0000000..fd2d3ad Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_P_err.png new file mode 100644 index 0000000..f581b69 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_P_flux.png new file mode 100644 index 0000000..5e4faa6 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_P_flux.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_SNRi.png new file mode 100644 index 0000000..da48cb3 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_SNRp.png new file mode 100644 index 0000000..2b5dd04 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_wae_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_crop_region.png b/plots/NGC1068_x274020/NGC1068_FOC_crop_region.png index 214919b..d05c2cb 100755 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_crop_region.png and b/plots/NGC1068_x274020/NGC1068_FOC_crop_region.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 73ea314..6e70762 100755 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 5d961f6..97b246c 100755 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/src/FOC_reduction.py b/src/FOC_reduction.py index 44a27ac..f095195 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -17,11 +17,11 @@ from lib.convex_hull import image_hull 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', @@ -60,9 +60,9 @@ def main(): # 'x3995202r_c0f.fits','x3995206r_c0f.fits'] # globals()['plots_folder'] = "../plots/PG1630+377_x39510/" - globals()['data_folder'] = "../data/IC5063_x3nl030/" - infiles = ['x3nl0301r_c0f.fits','x3nl0302r_c0f.fits','x3nl0303r_c0f.fits'] - globals()['plots_folder'] = "../plots/IC5063_x3nl030/" +# globals()['data_folder'] = "../data/IC5063_x3nl030/" +# infiles = ['x3nl0301r_c0f.fits','x3nl0302r_c0f.fits','x3nl0303r_c0f.fits'] +# globals()['plots_folder'] = "../plots/IC5063_x3nl030/" # globals()['data_folder'] = "../data/MKN3_x3nl010/" # infiles = ['x3nl0101r_c0f.fits','x3nl0102r_c0f.fits','x3nl0103r_c0f.fits'] @@ -77,10 +77,6 @@ def main(): # infiles = ['x3nl0201r_c0f.fits','x3nl0202r_c0f.fits','x3nl0203r_c0f.fits'] # globals()['plots_folder'] = "../plots/MKN78_x3nl020/" -# globals()['data_folder'] = "../data/PictorA_x25d040/" -# infiles = ['x25d0401t_c0f.fits','x25d0402t_c0f.fits','x25d0403t_c0f.fits'] -# globals()['plots_folder'] = "../plots/PictorA_x25d040/" - # globals()['data_folder'] = "../data/3C273_x0u20/" # infiles = ['x0u20101t_c0f.fits','x0u20102t_c0f.fits','x0u20103t_c0f.fits','x0u20104t_c0f.fits','x0u20105t_c0f.fits','x0u20106t_c0f.fits','x0u20201t_c0f.fits','x0u20202t_c0f.fits','x0u20203t_c0f.fits','x0u20204t_c0f.fits','x0u20205t_c0f.fits','x0u20206t_c0f.fits','x0u20301t_c0f.fits','x0u20302t_c0f.fits','x0u20303t_c0f.fits','x0u20304t_c0f.fits','x0u20305t_c0f.fits','x0u20306t_c0f.fits'] # globals()['plots_folder'] = "../plots/3C273_x0u20/" @@ -116,11 +112,11 @@ def main(): rotate_stokes = True #rotation to North convention can give erroneous results rotate_data = False #rotation to North convention can give erroneous results # Polarization map output - figname = 'IC5063_FOC' #target/intrument name - figtype = '_combine_FWHM020_pol' #additionnal informations + figname = 'NGC1068_FOC' #target/intrument name + figtype = '_combine_FWHM020_wae' #additionnal informations SNRp_cut = 3. #P measurments with SNR>3 - SNRi_cut = 70. #I measurments with SNR>30, which implies an uncertainty in P of 4.7%. - step_vec = 0 #plot all vectors in the array. if step_vec = 2, then every other vector will be plotted + 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 # if step_vec = 0 then all vectors are displayed at full length ##### Pipeline start @@ -163,7 +159,7 @@ def main(): else: vertex = np.array([0.,0.,data_array.shape[2],data_array.shape[2]]) shape = np.array([vertex[1]-vertex[0],vertex[3]-vertex[2]]) - rectangle = [vertex[2], vertex[0], shape[1], shape[0], 0., 'w'] + rectangle = [vertex[2], vertex[0], shape[1], shape[0], 0., 'g'] # Rotate data to have North up ref_header = deepcopy(headers[0]) @@ -186,7 +182,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 = 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, dP_dtheta, dPA_dtheta = 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 @@ -199,7 +195,7 @@ def main(): rectangle[4] = alpha 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, -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, 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, dP_dtheta, dPA_dtheta, headers) ## Step 4: # crop to desired region of interest (roi) diff --git a/src/lib/plots.py b/src/lib/plots.py index 430e274..cd918d5 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -93,16 +93,16 @@ def plot_obs(data_array, headers, shape=None, vmin=0., vmax=6., rectangle=None, exptime = headers[i]['exptime'] filt = headers[i]['filtnam1'] #plots - im = ax.imshow(data, vmin=vmin, vmax=vmax, origin='lower') + im = ax.imshow(data, vmin=vmin, vmax=vmax, origin='lower', cmap='gray') if not(rectangle is None): x, y, width, height, angle, color = rectangle[i] ax.add_patch(Rectangle((x, y), width, height, angle=angle, edgecolor=color, fill=False)) #position of centroid - ax.plot([data.shape[1]/2, data.shape[1]/2], [0,data.shape[0]-1], lw=1, - color='black') - ax.plot([0,data.shape[1]-1], [data.shape[1]/2, data.shape[1]/2], lw=1, - color='black') + ax.plot([data.shape[1]/2, data.shape[1]/2], [0,data.shape[0]-1], '--', lw=1, + color='grey', alpha=0.5) + ax.plot([0,data.shape[1]-1], [data.shape[1]/2, data.shape[1]/2], '--', lw=1, + color='grey', alpha=0.5) ax.annotate(instr+":"+rootname,color='white',fontsize=5,xy=(0.02, 0.95), xycoords='axes fraction') ax.annotate(filt,color='white',fontsize=10,xy=(0.02, 0.02), @@ -110,12 +110,12 @@ def plot_obs(data_array, headers, shape=None, vmin=0., vmax=6., rectangle=None, ax.annotate(exptime,color='white',fontsize=5,xy=(0.80, 0.02), xycoords='axes fraction') - fig.subplots_adjust(hspace=0, wspace=0, right=0.85) + fig.subplots_adjust(hspace=0.01, wspace=0.01, right=0.85) cbar_ax = fig.add_axes([0.9, 0.12, 0.02, 0.75]) - fig.colorbar(im, cax=cbar_ax) + fig.colorbar(im, cax=cbar_ax, label=r'$Counts \cdot s^{-1}$') if not (savename is None): - fig.suptitle(savename) + #fig.suptitle(savename) fig.savefig(plots_folder+savename+".png",bbox_inches='tight') plt.show() return 0 @@ -149,121 +149,27 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""): fig = plt.figure(figsize=(30,10)) ax = fig.add_subplot(131, projection=wcs) - im = ax.imshow(stkI, origin='lower') + im = ax.imshow(stkI, origin='lower', cmap='inferno') plt.colorbar(im) ax.set(xlabel="RA", ylabel="DEC", title=r"$I_{stokes}$") ax = fig.add_subplot(132, projection=wcs) - im = ax.imshow(stkQ, origin='lower') + im = ax.imshow(stkQ, origin='lower', cmap='inferno') plt.colorbar(im) ax.set(xlabel="RA", ylabel="DEC", title=r"$Q_{stokes}$") ax = fig.add_subplot(133, projection=wcs) - im = ax.imshow(stkU, origin='lower') + im = ax.imshow(stkU, origin='lower', cmap='inferno') plt.colorbar(im) ax.set(xlabel="RA", ylabel="DEC", title=r"$U_{stokes}$") if not (savename is None): - fig.suptitle(savename+"_IQU") + #fig.suptitle(savename+"_IQU") fig.savefig(plots_folder+savename+"_IQU.png",bbox_inches='tight') plt.show() return 0 -class crop_map(object): - """ - Class to interactively crop I_stokes map to desired Region of Interest - """ - def __init__(self, Stokes, data_mask, SNRp_cut=3., SNRi_cut=30.): - #Get data - stkI = Stokes[np.argmax([Stokes[i].header['datatype']=='I_stokes' for i in range(len(Stokes))])] - stk_cov = Stokes[np.argmax([Stokes[i].header['datatype']=='IQU_cov_matrix' for i in range(len(Stokes))])] - pol = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_deg_debiased' for i in range(len(Stokes))])] - pol_err = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_deg_err' for i in range(len(Stokes))])] - self.Stokes = Stokes - self.data_mask = data_mask - - wcs = WCS(Stokes[0]).deepcopy() - - #Compute SNR and apply cuts - pol.data[pol.data == 0.] = np.nan - SNRp = pol.data/pol_err.data - SNRp[np.isnan(SNRp)] = 0. - pol.data[SNRp < SNRp_cut] = np.nan - SNRi = stkI.data/np.sqrt(stk_cov.data[0,0]) - SNRi[np.isnan(SNRi)] = 0. - pol.data[SNRi < SNRi_cut] = np.nan - - convert_flux = Stokes[0].header['photflam'] - - #Plot the map - plt.rcParams.update({'font.size': 16}) - self.fig = plt.figure(figsize=(15,15)) - self.ax = self.fig.add_subplot(111, projection=wcs) - self.ax.set_facecolor('k') - self.fig.subplots_adjust(hspace=0, wspace=0, right=0.9) - cbar_ax = self.fig.add_axes([0.95, 0.12, 0.01, 0.75]) - - self.extent = [-stkI.data.shape[1]/2.,stkI.data.shape[1]/2.,-stkI.data.shape[0]/2.,stkI.data.shape[0]/2.] - self.center = [stkI.data.shape[1]/2.,stkI.data.shape[0]/2.] - - vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux) - im = self.ax.imshow(stkI.data*convert_flux,extent=self.extent, 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 = self.ax.contour(SNRi, extent=self.extent, levels=levelsI, colors='grey', linewidths=0.5) - - fontprops = fm.FontProperties(size=16) - px_size = wcs.wcs.get_cdelt()[0] - 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) - - self.ax.set_title( - "Click and drag to crop to desired Region of Interest.\n" - "Press 'c' to toggle the selector on and off") - - def onselect_crop(self, eclick, erelease) -> None: - # Obtain (xmin, xmax, ymin, ymax) values - self.RSextent = self.rect_selector.extents - self.RScenter = [self.center[i]+self.rect_selector.center[i] for i in range(2)] - - # Zoom to selection - print("CROP TO : ",self.RSextent) - print("CENTER : ",self.RScenter) - #self.ax.set_xlim(self.extent[0], self.extent[1]) - #self.ax.set_ylim(self.extent[2], self.extent[3]) - - def run(self) -> None: - self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, - drawtype='box', button=[1], interactive=True) - #self.fig.canvas.mpl_connect('key_press_event', self.toggle_selector) - plt.show() - - def crop(self): - Stokes_crop = copy.deepcopy(self.Stokes) - # Data sets to crop - stkI = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='I_stokes' for i in range(len(Stokes_crop))])] - stkQ = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Q_stokes' for i in range(len(Stokes_crop))])] - stkU = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='U_stokes' for i in range(len(Stokes_crop))])] - stk_cov = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='IQU_cov_matrix' for i in range(len(Stokes_crop))])] - pol = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_deg_debiased' for i in range(len(Stokes_crop))])] - pol_err = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_deg_err' for i in range(len(Stokes_crop))])] - pang = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_ang' for i in range(len(Stokes_crop))])] - pang_err = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_ang_err' for i in range(len(Stokes_crop))])] - # Crop datasets - vertex = [int(self.RScenter[0]+self.RSextent[0]), int(self.RScenter[0]+self.RSextent[1]), int(self.RScenter[1]+self.RSextent[2]), int(self.RScenter[1]+self.RSextent[3])] - for dataset in [stkI, stkQ, stkU, pol, pol_err, pang, pang_err]: - dataset.data = dataset.data[vertex[2]:vertex[3], vertex[0]:vertex[1]] - data_mask = self.data_mask[vertex[2]:vertex[3], vertex[0]:vertex[1]] - new_stk_cov = np.zeros((3,3,stkI.data.shape[0],stkI.data.shape[1])) - for i in range(3): - for j in range(3): - new_stk_cov[i,j] = stk_cov.data[i,j][vertex[2]:vertex[3], vertex[0]:vertex[1]] - stk_cov.data = new_stk_cov - - return Stokes_crop, data_mask - - def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_cut=30., step_vec=1, savename=None, plots_folder="", display=None): """ @@ -366,7 +272,9 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c im = ax.imshow(stkI.data*convert_flux, 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) + print("SNRi contour levels : ", levelsI) cont = ax.contour(SNRi, levels=levelsI, colors='grey', linewidths=0.5) + #ax.clabel(cont,inline=True,fontsize=6) elif display.lower() in ['pol_flux']: # Display polarisation flux pf_mask = (stkI.data > 0.) * (pol.data > 0.) @@ -673,4 +581,99 @@ class overplot_maps(align_maps): self.align() if self.aligned: self.overplot(other_levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename) - plt.show(block=True) \ No newline at end of file + plt.show(block=True) + + +class crop_map(object): + """ + Class to interactively crop I_stokes map to desired Region of Interest + """ + def __init__(self, Stokes, data_mask, SNRp_cut=3., SNRi_cut=30.): + #Get data + stkI = Stokes[np.argmax([Stokes[i].header['datatype']=='I_stokes' for i in range(len(Stokes))])] + stk_cov = Stokes[np.argmax([Stokes[i].header['datatype']=='IQU_cov_matrix' for i in range(len(Stokes))])] + pol = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_deg_debiased' for i in range(len(Stokes))])] + pol_err = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_deg_err' for i in range(len(Stokes))])] + self.Stokes = Stokes + self.data_mask = data_mask + + wcs = WCS(Stokes[0]).deepcopy() + + #Compute SNR and apply cuts + pol.data[pol.data == 0.] = np.nan + SNRp = pol.data/pol_err.data + SNRp[np.isnan(SNRp)] = 0. + pol.data[SNRp < SNRp_cut] = np.nan + SNRi = stkI.data/np.sqrt(stk_cov.data[0,0]) + SNRi[np.isnan(SNRi)] = 0. + pol.data[SNRi < SNRi_cut] = np.nan + + convert_flux = Stokes[0].header['photflam'] + + #Plot the map + plt.rcParams.update({'font.size': 16}) + self.fig = plt.figure(figsize=(15,15)) + self.ax = self.fig.add_subplot(111, projection=wcs) + self.ax.set_facecolor('k') + self.fig.subplots_adjust(hspace=0, wspace=0, right=0.9) + cbar_ax = self.fig.add_axes([0.95, 0.12, 0.01, 0.75]) + + self.extent = [-stkI.data.shape[1]/2.,stkI.data.shape[1]/2.,-stkI.data.shape[0]/2.,stkI.data.shape[0]/2.] + self.center = [stkI.data.shape[1]/2.,stkI.data.shape[0]/2.] + + vmin, vmax = 0., np.max(stkI.data[stkI.data > 0.]*convert_flux) + im = self.ax.imshow(stkI.data*convert_flux,extent=self.extent, 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 = self.ax.contour(SNRi, extent=self.extent, levels=levelsI, colors='grey', linewidths=0.5) + + fontprops = fm.FontProperties(size=16) + px_size = wcs.wcs.get_cdelt()[0] + 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) + + self.ax.set_title( + "Click and drag to crop to desired Region of Interest.\n" + "Press 'c' to toggle the selector on and off") + + def onselect_crop(self, eclick, erelease) -> None: + # Obtain (xmin, xmax, ymin, ymax) values + self.RSextent = self.rect_selector.extents + self.RScenter = [self.center[i]+self.rect_selector.center[i] for i in range(2)] + + # Zoom to selection + print("CROP TO : ",self.RSextent) + print("CENTER : ",self.RScenter) + #self.ax.set_xlim(self.extent[0], self.extent[1]) + #self.ax.set_ylim(self.extent[2], self.extent[3]) + + def run(self) -> None: + self.rect_selector = RectangleSelector(self.ax, self.onselect_crop, + drawtype='box', button=[1], interactive=True) + #self.fig.canvas.mpl_connect('key_press_event', self.toggle_selector) + plt.show() + + def crop(self): + Stokes_crop = copy.deepcopy(self.Stokes) + # Data sets to crop + stkI = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='I_stokes' for i in range(len(Stokes_crop))])] + stkQ = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Q_stokes' for i in range(len(Stokes_crop))])] + stkU = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='U_stokes' for i in range(len(Stokes_crop))])] + stk_cov = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='IQU_cov_matrix' for i in range(len(Stokes_crop))])] + pol = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_deg_debiased' for i in range(len(Stokes_crop))])] + pol_err = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_deg_err' for i in range(len(Stokes_crop))])] + pang = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_ang' for i in range(len(Stokes_crop))])] + pang_err = Stokes_crop[np.argmax([Stokes_crop[i].header['datatype']=='Pol_ang_err' for i in range(len(Stokes_crop))])] + # Crop datasets + vertex = [int(self.RScenter[0]+self.RSextent[0]), int(self.RScenter[0]+self.RSextent[1]), int(self.RScenter[1]+self.RSextent[2]), int(self.RScenter[1]+self.RSextent[3])] + for dataset in [stkI, stkQ, stkU, pol, pol_err, pang, pang_err]: + dataset.data = dataset.data[vertex[2]:vertex[3], vertex[0]:vertex[1]] + data_mask = self.data_mask[vertex[2]:vertex[3], vertex[0]:vertex[1]] + new_stk_cov = np.zeros((3,3,stkI.data.shape[0],stkI.data.shape[1])) + for i in range(3): + for j in range(3): + new_stk_cov[i,j] = stk_cov.data[i,j][vertex[2]:vertex[3], vertex[0]:vertex[1]] + stk_cov.data = new_stk_cov + + return Stokes_crop, data_mask + diff --git a/src/lib/reduction.py b/src/lib/reduction.py index cfb65ef..d745e2a 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -246,21 +246,21 @@ def crop_array(data_array, headers, error_array=None, step=5, null_val=None, new_shape = np.array([v_array[1]-v_array[0],v_array[3]-v_array[2]]) rectangle = [v_array[2], v_array[0], new_shape[1], new_shape[0], 0., 'b'] if display: - fig, ax = plt.subplots() + fig, ax = plt.subplots(figsize=(10,10)) data = data_array[0] instr = headers[0]['instrume'] rootname = headers[0]['rootname'] exptime = headers[0]['exptime'] filt = headers[0]['filtnam1'] #plots - im = ax.imshow(data, vmin=data.min(), vmax=data.max(), origin='lower') + im = ax.imshow(data, vmin=data.min(), vmax=data.max(), origin='lower', cmap='gray') x, y, width, height, angle, color = rectangle ax.add_patch(Rectangle((x, y),width,height,edgecolor=color,fill=False)) #position of centroid - ax.plot([data.shape[1]/2, data.shape[1]/2], [0,data.shape[0]-1], lw=1, - color='black') - ax.plot([0,data.shape[1]-1], [data.shape[1]/2, data.shape[1]/2], lw=1, - color='black') + ax.plot([data.shape[1]/2, data.shape[1]/2], [0,data.shape[0]-1], '--', lw=1, + color='grey', alpha=0.3) + ax.plot([0,data.shape[1]-1], [data.shape[1]/2, data.shape[1]/2], '--', lw=1, + color='grey', alpha=0.3) ax.annotate(instr+":"+rootname, color='white', fontsize=5, xy=(0.02, 0.95), xycoords='axes fraction') ax.annotate(filt, color='white', fontsize=10, xy=(0.02, 0.02), @@ -273,10 +273,10 @@ def crop_array(data_array, headers, error_array=None, step=5, null_val=None, fig.subplots_adjust(hspace=0, wspace=0, right=0.85) cbar_ax = fig.add_axes([0.9, 0.12, 0.02, 0.75]) - fig.colorbar(im, cax=cbar_ax) + fig.colorbar(im, cax=cbar_ax, label=r'$Counts \cdot s^{-1}$') if not(savename is None): - fig.suptitle(savename+'_'+filt+'_crop_region') + #fig.suptitle(savename+'_'+filt+'_crop_region') fig.savefig(plots_folder+savename+'_'+filt+'_crop_region.png', bbox_inches='tight') plot_obs(data_array, headers, vmin=data_array.min(), @@ -486,7 +486,7 @@ def get_error(data_array, headers, sub_shape=(15,15), display=False, plt.legend() if not(savename is None): - fig.suptitle(savename+"_background_flux") + #fig.suptitle(savename+"_background_flux") fig.savefig(plots_folder+savename+"_background_flux.png", bbox_inches='tight') vmin = np.min(np.log10(data[data > 0.])) @@ -584,7 +584,11 @@ def rebin_array(data_array, error_array, headers, pxsize, scale, # Propagate error rms_image = np.sqrt(bin_ndarray(image**2, new_shape=new_shape, operation='average')) - new_error = np.sqrt(Dxy[0]*Dxy[1])*bin_ndarray(error, + if operation.lower() in ["mean", "average", "avg"]: + new_error = 1./np.sqrt(Dxy[0]*Dxy[1])*bin_ndarray(error, + new_shape=new_shape, operation='average') + else: + new_error = np.sqrt(Dxy[0]*Dxy[1])*bin_ndarray(error, new_shape=new_shape, operation='average') rebinned_error.append(np.sqrt(rms_image**2 + new_error**2)) @@ -1083,11 +1087,12 @@ def compute_Stokes(data_array, error_array, data_mask, headers, pol_eff[1] = pol_efficiency['pol60'] pol_eff[2] = pol_efficiency['pol120'] - # Orientation and error for each polarizer ## THIS IS WHERE WE IMPLEMENT THE ERROR THAT IS GOING WRONG + # Orientation and error for each polarizer # POL0 = 0deg, POL60 = 60deg, POL120=120deg theta = np.array([180.*np.pi/180., 60.*np.pi/180., 120.*np.pi/180.]) # Uncertainties on the orientation of the polarizers' axes taken to be 3deg (see Nota et. al 1996, p36; Robinson & Thomson 1995) sigma_theta = np.array([3.*np.pi/180., 3.*np.pi/180., 3.*np.pi/180.]) + fmax = np.finfo(np.float64).max pol_flux = 2.*np.array([pol0/transmit[0], pol60/transmit[1], pol120/transmit[2]]) pol_cov = np.array([[4.*pol_cov_m[i,j]/(transmit[i]*transmit[j]) for j in range(pol_cov_m.shape[1])] for i in range(pol_cov_m.shape[0])]) @@ -1130,19 +1135,42 @@ def compute_Stokes(data_array, error_array, data_mask, headers, dI_dtheta1 = 2.*pol_eff[0]/A*(pol_eff[2]*np.cos(-2.*theta[2]+2.*theta[0])*(pol_flux[1]-I_stokes) - pol_eff[1]*np.cos(-2.*theta[0]+2.*theta[1])*(pol_flux[2]-I_stokes)) 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 derivative of each polarization component with respect to the polarizer orientation + mask = I_stokes > 0. + dP_dtheta1 = np.ones(I_stokes.shape)*fmax + dP_dtheta2 = np.ones(I_stokes.shape)*fmax + dP_dtheta3 = np.ones(I_stokes.shape)*fmax + dP_dtheta1[mask] = 1./(I_stokes[mask]*np.sqrt(Q_stokes[mask]**2 + U_stokes[mask]**2))*(Q_stokes[mask]*dQ_dtheta1[mask] + U_stokes[mask]*dQ_dtheta1[mask]) - 1./I_stokes[mask]**2 * dI_dtheta1[mask] + dP_dtheta2[mask] = 1./(I_stokes[mask]*np.sqrt(Q_stokes[mask]**2 + U_stokes[mask]**2))*(Q_stokes[mask]*dQ_dtheta2[mask] + U_stokes[mask]*dQ_dtheta2[mask]) - 1./I_stokes[mask]**2 * dI_dtheta2[mask] + dP_dtheta3[mask] = 1./(I_stokes[mask]*np.sqrt(Q_stokes[mask]**2 + U_stokes[mask]**2))*(Q_stokes[mask]*dQ_dtheta3[mask] + U_stokes[mask]*dQ_dtheta3[mask]) - 1./I_stokes[mask]**2 * dI_dtheta3[mask] + dP_dtheta = np.array([dP_dtheta1, dP_dtheta2, dP_dtheta3]) + + dPA_dtheta1 = np.ones(I_stokes.shape)*fmax + dPA_dtheta2 = np.ones(I_stokes.shape)*fmax + dPA_dtheta3 = np.ones(I_stokes.shape)*fmax + dPA_dtheta1[mask] = 90./np.pi*1./(U_stokes[mask]**2 + Q_stokes[mask]**2)*(Q_stokes[mask]*dU_dtheta1[mask] - U_stokes[mask]*dQ_dtheta1[mask]) + dPA_dtheta2[mask] = 90./np.pi*1./(U_stokes[mask]**2 + Q_stokes[mask]**2)*(Q_stokes[mask]*dU_dtheta2[mask] - U_stokes[mask]*dQ_dtheta2[mask]) + dPA_dtheta3[mask] = 90./np.pi*1./(U_stokes[mask]**2 + Q_stokes[mask]**2)*(Q_stokes[mask]*dU_dtheta3[mask] - U_stokes[mask]*dQ_dtheta3[mask]) + dPA_dtheta = np.array([dPA_dtheta1, dPA_dtheta2, dPA_dtheta3]) # Compute the uncertainty associated with the polarizers' orientation (see Kishimoto 1999) - s_I2_axis = (dI_dtheta1**2*sigma_theta[0]**2 + dI_dtheta2**2*sigma_theta[1]**2 + dI_dtheta3**2*sigma_theta[2]**2) - s_Q2_axis = (dQ_dtheta1**2*sigma_theta[0]**2 + dQ_dtheta2**2*sigma_theta[1]**2 + dQ_dtheta3**2*sigma_theta[2]**2) - s_U2_axis = (dU_dtheta1**2*sigma_theta[0]**2 + dU_dtheta2**2*sigma_theta[1]**2 + dU_dtheta3**2*sigma_theta[2]**2) + s_I2_axis = np.sum([dI_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))],axis=0) #(dI_dtheta1**2*sigma_theta[0]**2 + dI_dtheta2**2*sigma_theta[1]**2 + dI_dtheta3**2*sigma_theta[2]**2) + s_Q2_axis = np.sum([dQ_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))],axis=0) #(dQ_dtheta1**2*sigma_theta[0]**2 + dQ_dtheta2**2*sigma_theta[1]**2 + dQ_dtheta3**2*sigma_theta[2]**2) + s_U2_axis = np.sum([dU_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))],axis=0) #(dU_dtheta1**2*sigma_theta[0]**2 + dU_dtheta2**2*sigma_theta[1]**2 + dU_dtheta3**2*sigma_theta[2]**2) - # Add quadratically the uncertainty to the Stokes covariance matrix ## THIS IS WHERE THE PROBLEMATIC UNCERTAINTY IS ADDED TO THE PIPELINE + # Add quadratically the uncertainty to the Stokes covariance matrix Stokes_cov[0,0] += s_I2_axis Stokes_cov[1,1] += s_Q2_axis Stokes_cov[2,2] += s_U2_axis @@ -1158,10 +1186,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 + return I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta -def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers): +def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta, headers): """ Compute the polarization degree (in %) and angle (in deg) and their respective errors from given Stokes parameters. @@ -1215,13 +1243,24 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers): print("WARNING : found {0:d} pixels for which P > 1".format(P[P>1.].size)) #Associated errors + sigma_theta = np.array([3.*np.pi/180., 3.*np.pi/180., 3.*np.pi/180.]) fmax = np.finfo(np.float64).max - s_P = np.ones(I_stokes.shape)*fmax - s_P[mask] = (1/I_stokes[mask])*np.sqrt((Q_stokes[mask]**2*Stokes_cov[1,1][mask] + U_stokes[mask]**2*Stokes_cov[2,2][mask] + 2.*Q_stokes[mask]*U_stokes[mask]*Stokes_cov[1,2][mask])/(Q_stokes[mask]**2 + U_stokes[mask]**2) + ((Q_stokes[mask]/I_stokes[mask])**2 + (U_stokes[mask]/I_stokes[mask])**2)*Stokes_cov[0,0][mask] - 2.*(Q_stokes[mask]/I_stokes[mask])*Stokes_cov[0,1][mask] - 2.*(U_stokes[mask]/I_stokes[mask])*Stokes_cov[0,2][mask]) - s_P[np.isnan(s_P)] = fmax s_PA = np.ones(I_stokes.shape)*fmax + + # Propagate previously computed errors + s_P[mask] = (1/I_stokes[mask])*np.sqrt((Q_stokes[mask]**2*Stokes_cov[1,1][mask] + U_stokes[mask]**2*Stokes_cov[2,2][mask] + 2.*Q_stokes[mask]*U_stokes[mask]*Stokes_cov[1,2][mask])/(Q_stokes[mask]**2 + U_stokes[mask]**2) + ((Q_stokes[mask]/I_stokes[mask])**2 + (U_stokes[mask]/I_stokes[mask])**2)*Stokes_cov[0,0][mask] - 2.*(Q_stokes[mask]/I_stokes[mask])*Stokes_cov[0,1][mask] - 2.*(U_stokes[mask]/I_stokes[mask])*Stokes_cov[0,2][mask]) s_PA[mask] = (90./(np.pi*(Q_stokes[mask]**2 + U_stokes[mask]**2)))*np.sqrt(U_stokes[mask]**2*Stokes_cov[1,1][mask] + Q_stokes[mask]**2*Stokes_cov[2,2][mask] - 2.*Q_stokes[mask]*U_stokes[mask]*Stokes_cov[1,2][mask]) + # Compute the uncertainty associated with the polarizers' orientation (see Kishimoto 1999) + s_P_axis = np.sqrt(np.sum([dP_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))], axis=0)) + s_PA_axis = np.sqrt(np.sum([dPA_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))], axis=0)) + # Add axis uncertainties quadratically +# with warnings.catch_warnings(record=True) as w: +# s_P[mask] = np.sqrt(s_P[mask]**2 + s_P_axis[mask]**2) +# s_PA[mask] = np.sqrt(s_PA[mask]**2 + s_PA_axis[mask]**2) + + s_P[np.isnan(s_P)] = fmax + s_PA[np.isnan(s_PA)] = fmax # Catch expected "OverflowWarning" as wrong pixel have an overflowing error with warnings.catch_warnings(record=True) as w: diff --git a/src/lib/test_overplot.py b/src/lib/test_overplot.py index 8aead6c..1e04e69 100755 --- a/src/lib/test_overplot.py +++ b/src/lib/test_overplot.py @@ -12,9 +12,9 @@ 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_maps(Stokes_UV, Stokes_18GHz) -A.plot(levels=levels18GHz, SNRp_cut=6.0, SNRi_cut=180.0, savename='../../plots/IC5063_x3nl030/18GHz_overplot.png') +A.plot(levels=levels18GHz, SNRp_cut=3.0, SNRi_cut=50.0, savename='../../plots/IC5063_x3nl030/18GHz_overplot.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_maps(Stokes_UV, Stokes_24GHz) -B.plot(levels=levels24GHz, SNRp_cut=6.0, SNRi_cut=180.0, savename='../../plots/IC5063_x3nl030/24GHz_overplot.png') +B.plot(levels=levels24GHz, SNRp_cut=3.0, SNRi_cut=50.0, savename='../../plots/IC5063_x3nl030/24GHz_overplot.png')