Tibeuleu/FOC_Reduction
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

polarisation -> polarization everywhere

Browse Source
This commit is contained in:
Thibault Barnouin
2024-06-04 11:59:33 +02:00
parent 5e27f36869
commit 97479f90fb
5 changed files with 90 additions and 96 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
package/FOC_reduction.py
View File

@@ -165,7 +165,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, SNRi_cut=None) I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, SNRi_cut=None)
I_bkg, Q_bkg, U_bkg, S_cov_bkg, _, _ = proj_red.rotate_Stokes(I_bkg, Q_bkg, U_bkg, S_cov_bkg, np.array(True).reshape(1, 1), headers, SNRi_cut=None) I_bkg, Q_bkg, U_bkg, S_cov_bkg, _, _ = proj_red.rotate_Stokes(I_bkg, Q_bkg, U_bkg, S_cov_bkg, np.array(True).reshape(1, 1), headers, SNRi_cut=None)
# Compute polarimetric parameters (polarisation degree and angle). # 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, headers)
P_bkg, debiased_P_bkg, s_P_bkg, s_P_P_bkg, PA_bkg, s_PA_bkg, s_PA_P_bkg = proj_red.compute_pol(I_bkg, Q_bkg, U_bkg, S_cov_bkg, headers) P_bkg, debiased_P_bkg, s_P_bkg, s_P_P_bkg, PA_bkg, s_PA_bkg, s_PA_P_bkg = proj_red.compute_pol(I_bkg, Q_bkg, U_bkg, S_cov_bkg, headers)
@@ -194,28 +194,28 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
I_bkg[0, 0]*headers[0]['photflam'], np.sqrt(S_cov_bkg[0, 0][0, 0])*headers[0]['photflam'], 2, out=int))) I_bkg[0, 0]*headers[0]['photflam'], np.sqrt(S_cov_bkg[0, 0][0, 0])*headers[0]['photflam'], 2, out=int)))
print("P_bkg = {0:.1f} ± {1:.1f} %".format(debiased_P_bkg[0, 0]*100., np.ceil(s_P_bkg[0, 0]*1000.)/10.)) print("P_bkg = {0:.1f} ± {1:.1f} %".format(debiased_P_bkg[0, 0]*100., np.ceil(s_P_bkg[0, 0]*1000.)/10.))
print("PA_bkg = {0:.1f} ± {1:.1f} °".format(princ_angle(PA_bkg[0, 0]), princ_angle(np.ceil(s_PA_bkg[0, 0]*10.)/10.))) print("PA_bkg = {0:.1f} ± {1:.1f} °".format(princ_angle(PA_bkg[0, 0]), princ_angle(np.ceil(s_PA_bkg[0, 0]*10.)/10.)))
# Plot polarisation map (Background is either total Flux, Polarization degree or Polarization degree error). # Plot polarization map (Background is either total Flux, Polarization degree or Polarization degree error).
if px_scale.lower() not in ['full', 'integrate'] and not interactive: if px_scale.lower() not in ['full', 'integrate'] and not interactive:
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim,
step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname]), plots_folder=plots_folder) step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname]), plots_folder=plots_folder)
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec,
vec_scale=vec_scale, savename="_".join([figname, "I"]), plots_folder=plots_folder, display='Intensity') vec_scale=vec_scale, savename="_".join([figname, "I"]), plots_folder=plots_folder, display='Intensity')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec,
vec_scale=vec_scale, savename="_".join([figname, "P_flux"]), plots_folder=plots_folder, display='Pol_Flux') vec_scale=vec_scale, savename="_".join([figname, "P_flux"]), plots_folder=plots_folder, display='Pol_Flux')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec,
vec_scale=vec_scale, savename="_".join([figname, "P"]), plots_folder=plots_folder, display='Pol_deg') vec_scale=vec_scale, savename="_".join([figname, "P"]), plots_folder=plots_folder, display='Pol_deg')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec,
vec_scale=vec_scale, savename="_".join([figname, "PA"]), plots_folder=plots_folder, display='Pol_ang') vec_scale=vec_scale, savename="_".join([figname, "PA"]), plots_folder=plots_folder, display='Pol_ang')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec,
vec_scale=vec_scale, savename="_".join([figname, "I_err"]), plots_folder=plots_folder, display='I_err') vec_scale=vec_scale, savename="_".join([figname, "I_err"]), plots_folder=plots_folder, display='I_err')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec,
vec_scale=vec_scale, savename="_".join([figname, "P_err"]), plots_folder=plots_folder, display='Pol_deg_err') vec_scale=vec_scale, savename="_".join([figname, "P_err"]), plots_folder=plots_folder, display='Pol_deg_err')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec,
vec_scale=vec_scale, savename="_".join([figname, "SNRi"]), plots_folder=plots_folder, display='SNRi') vec_scale=vec_scale, savename="_".join([figname, "SNRi"]), plots_folder=plots_folder, display='SNRi')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec,
vec_scale=vec_scale, savename="_".join([figname, "SNRp"]), plots_folder=plots_folder, display='SNRp') vec_scale=vec_scale, savename="_".join([figname, "SNRp"]), plots_folder=plots_folder, display='SNRp')
elif not interactive: elif not interactive:
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut,
savename=figname, plots_folder=plots_folder, display='integrate') savename=figname, plots_folder=plots_folder, display='integrate')
elif px_scale.lower() not in ['full', 'integrate']: elif px_scale.lower() not in ['full', 'integrate']:
proj_plots.pol_map(Stokes_test, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim) proj_plots.pol_map(Stokes_test, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim)

8
package/lib/fits.py
View File

@@ -144,10 +144,10 @@ def save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P,
header['targname'] = (ref_header['targname'], 'Target name') header['targname'] = (ref_header['targname'], 'Target name')
header['orientat'] = (ref_header['orientat'], 'Angle between North and the y-axis of the image') header['orientat'] = (ref_header['orientat'], 'Angle between North and the y-axis of the image')
header['filename'] = (filename, 'Original filename') header['filename'] = (filename, 'Original filename')
header['P_int'] = (ref_header['P_int'], 'Integrated polarisation degree') header['P_int'] = (ref_header['P_int'], 'Integrated polarization degree')
header['P_int_err'] = (ref_header['P_int_err'], 'Integrated polarisation degree error') header['P_int_err'] = (ref_header['P_int_err'], 'Integrated polarization degree error')
header['PA_int'] = (ref_header['PA_int'], 'Integrated polarisation angle') header['PA_int'] = (ref_header['PA_int'], 'Integrated polarization angle')
header['PA_int_err'] = (ref_header['PA_int_err'], 'Integrated polarisation angle error') header['PA_int_err'] = (ref_header['PA_int_err'], 'Integrated polarization angle error')
# Crop Data to mask # Crop Data to mask
if data_mask.shape != (1, 1): if data_mask.shape != (1, 1):

82
package/lib/plots.py
View File

@@ -9,8 +9,8 @@ prototypes :
- plot_Stokes(Stokes, savename, plots_folder) - plot_Stokes(Stokes, savename, plots_folder)
Plot the I/Q/U maps from the Stokes HDUList. Plot the I/Q/U maps from the Stokes HDUList.
- polarisation_map(Stokes, data_mask, rectangle, SNRp_cut, SNRi_cut, step_vec, savename, plots_folder, display) -> fig, ax - polarization_map(Stokes, data_mask, rectangle, SNRp_cut, SNRi_cut, step_vec, savename, plots_folder, display) -> fig, ax
Plots polarisation map of polarimetric parameters saved in an HDUList. Plots polarization map of polarimetric parameters saved in an HDUList.
class align_maps(map, other_map, **kwargs) class align_maps(map, other_map, **kwargs)
Class to interactively align maps with different WCS. Class to interactively align maps with different WCS.
@@ -22,13 +22,13 @@ prototypes :
Class inherited from align_maps to overplot chandra data as contours. Class inherited from align_maps to overplot chandra data as contours.
class overplot_pol(align_maps) class overplot_pol(align_maps)
Class inherited from align_maps to overplot UV polarisation vectors on other maps. Class inherited from align_maps to overplot UV polarization vectors on other maps.
class crop_map(hdul, fig, ax) class crop_map(hdul, fig, ax)
Class to interactively crop a region of interest of a HDUList. Class to interactively crop a region of interest of a HDUList.
class crop_Stokes(crop_map) class crop_Stokes(crop_map)
Class inherited from crop_map to work on polarisation maps. Class inherited from crop_map to work on polarization maps.
class image_lasso_selector(img, fig, ax) class image_lasso_selector(img, fig, ax)
Class to interactively select part of a map to work on. Class to interactively select part of a map to work on.
@@ -37,7 +37,7 @@ prototypes :
Class to interactively simulate aperture integration. Class to interactively simulate aperture integration.
class pol_map(Stokes, SNRp_cut, SNRi_cut, selection) class pol_map(Stokes, SNRp_cut, SNRi_cut, selection)
Class to interactively study polarisation maps making use of the cropping and selecting tools. Class to interactively study polarization maps making use of the cropping and selecting tools.
""" """
from copy import deepcopy from copy import deepcopy
@@ -204,10 +204,10 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""):
return 0 return 0
def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_cut=3., def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_cut=3.,
flux_lim=None, step_vec=1, vec_scale=2., savename=None, plots_folder="", display="default"): flux_lim=None, step_vec=1, vec_scale=2., savename=None, plots_folder="", display="default"):
""" """
Plots polarisation map from Stokes HDUList. Plots polarization map from Stokes HDUList.
---------- ----------
Inputs: Inputs:
Stokes : astropy.io.fits.hdu.hdulist.HDUList Stokes : astropy.io.fits.hdu.hdulist.HDUList
@@ -230,12 +230,12 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
Defaults to None, limits are computed on the background value and the Defaults to None, limits are computed on the background value and the
maximum value in the cut. maximum value in the cut.
step_vec : int, optional step_vec : int, optional
Number of steps between each displayed polarisation vector. Number of steps between each displayed polarization vector.
If step_vec = 2, every other vector will be displayed. If step_vec = 2, every other vector will be displayed.
Defaults to 1 Defaults to 1
vec_scale : float, optional vec_scale : float, optional
Pixel length of displayed 100% polarisation vector. Pixel length of displayed 100% polarization vector.
If vec_scale = 2, a vector of 50% polarisation will be 1 pixel wide. If vec_scale = 2, a vector of 50% polarization will be 1 pixel wide.
Defaults to 2. Defaults to 2.
savename : str, optional savename : str, optional
Name of the figure the map should be saved to. If None, the map won't Name of the figure the map should be saved to. If None, the map won't
@@ -246,8 +246,8 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
be saved. Not used if savename is None. be saved. Not used if savename is None.
Defaults to current folder. Defaults to current folder.
display : str, optional display : str, optional
Choose the map to display between intensity (default), polarisation Choose the map to display between intensity (default), polarization
degree ('p', 'pol', 'pol_deg') or polarisation degree error ('s_p', degree ('p', 'pol', 'pol_deg') or polarization degree error ('s_p',
'pol_err', 'pol_deg_err'). 'pol_err', 'pol_deg_err').
Defaults to None (intensity). Defaults to None (intensity).
---------- ----------
@@ -295,12 +295,12 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
poldata[np.logical_not(mask)] = np.nan poldata[np.logical_not(mask)] = np.nan
pangdata[np.logical_not(mask)] = np.nan pangdata[np.logical_not(mask)] = np.nan
# Look for pixel of max polarisation # Look for pixel of max polarization
if np.isfinite(pol).any(): if np.isfinite(pol).any():
p_max = np.max(pol[np.isfinite(pol)]) p_max = np.max(pol[np.isfinite(pol)])
x_max, y_max = np.unravel_index(np.argmax(pol == p_max), pol.shape) x_max, y_max = np.unravel_index(np.argmax(pol == p_max), pol.shape)
else: else:
print("No pixel with polarisation information above requested SNR.") print("No pixel with polarization information above requested SNR.")
# Plot the map # Plot the map
plt.rcParams.update({'font.size': 10}) plt.rcParams.update({'font.size': 10})
@@ -325,7 +325,7 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
print("Total flux contour levels : ", levelsI) print("Total flux contour levels : ", levelsI)
ax.contour(stkI*convert_flux, levels=levelsI, colors='grey', linewidths=0.5) ax.contour(stkI*convert_flux, levels=levelsI, colors='grey', linewidths=0.5)
elif display.lower() in ['pol_flux']: elif display.lower() in ['pol_flux']:
# Display polarisation flux # Display polarization flux
display = 'pf' display = 'pf'
if flux_lim is None: if flux_lim is None:
if mask.sum() > 0.: if mask.sum() > 0.:
@@ -340,19 +340,19 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
print("Polarized flux contour levels : ", levelsPf) print("Polarized flux contour levels : ", levelsPf)
ax.contour(stkI*convert_flux*pol, levels=levelsPf, colors='grey', linewidths=0.5) ax.contour(stkI*convert_flux*pol, levels=levelsPf, colors='grey', linewidths=0.5)
elif display.lower() in ['p', 'pol', 'pol_deg']: elif display.lower() in ['p', 'pol', 'pol_deg']:
# Display polarisation degree map # Display polarization degree map
display = 'p' display = 'p'
vmin, vmax = 0., 100. vmin, vmax = 0., 100.
im = ax.imshow(pol*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) im = ax.imshow(pol*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.)
fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$P$ [%]") fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$P$ [%]")
elif display.lower() in ['pa', 'pang', 'pol_ang']: elif display.lower() in ['pa', 'pang', 'pol_ang']:
# Display polarisation degree map # Display polarization degree map
display = 'pa' display = 'pa'
vmin, vmax = 0., 180. vmin, vmax = 0., 180.
im = ax.imshow(princ_angle(pang), vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) im = ax.imshow(princ_angle(pang), vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.)
fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$\theta_P$ [°]") fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$\theta_P$ [°]")
elif display.lower() in ['s_p', 'pol_err', 'pol_deg_err']: elif display.lower() in ['s_p', 'pol_err', 'pol_deg_err']:
# Display polarisation degree error map # Display polarization degree error map
display = 's_p' display = 's_p'
if (SNRp > SNRp_cut).any(): if (SNRp > SNRp_cut).any():
vmin, vmax = 0., np.max([pol_err[SNRp > SNRp_cut].max(), 1.])*100. vmin, vmax = 0., np.max([pol_err[SNRp > SNRp_cut].max(), 1.])*100.
@@ -384,7 +384,7 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
im = ax.imshow(SNRi, aspect='equal', cmap='inferno', alpha=1.) im = ax.imshow(SNRi, aspect='equal', cmap='inferno', alpha=1.)
fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$I_{Stokes}/\sigma_{I}$") fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$I_{Stokes}/\sigma_{I}$")
elif display.lower() in ['snrp']: elif display.lower() in ['snrp']:
# Display polarisation degree signal-to-noise map # Display polarization degree signal-to-noise map
display = 'snrp' display = 'snrp'
vmin, vmax = 0., np.max(SNRp[np.isfinite(SNRp)]) vmin, vmax = 0., np.max(SNRp[np.isfinite(SNRp)])
if vmax*0.99 > SNRp_cut: if vmax*0.99 > SNRp_cut:
@@ -710,7 +710,7 @@ class overplot_radio(align_maps):
self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(10, 10), subplot_kw=dict(projection=self.wcs_UV)) self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(10, 10), subplot_kw=dict(projection=self.wcs_UV))
self.fig_overplot.subplots_adjust(hspace=0, wspace=0, bottom=0.1, left=0.1, top=0.8, right=1) self.fig_overplot.subplots_adjust(hspace=0, wspace=0, bottom=0.1, left=0.1, top=0.8, right=1)
# Display UV intensity map with polarisation vectors # Display UV intensity map with polarization vectors
vmin, vmax = stkI[np.isfinite(stkI)].max()/1e3*self.map_convert, stkI[np.isfinite(stkI)].max()*self.map_convert vmin, vmax = stkI[np.isfinite(stkI)].max()/1e3*self.map_convert, stkI[np.isfinite(stkI)].max()*self.map_convert
for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]:
try: try:
@@ -728,7 +728,7 @@ class overplot_radio(align_maps):
self.cbar = self.fig_overplot.colorbar(self.im, ax=self.ax_overplot, aspect=50, shrink=0.75, pad=0.025, self.cbar = self.fig_overplot.colorbar(self.im, ax=self.ax_overplot, aspect=50, shrink=0.75, pad=0.025,
label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit)) label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit))
# Display full size polarisation vectors # Display full size polarization vectors
if vec_scale is None: if vec_scale is None:
self.vec_scale = 2. self.vec_scale = 2.
pol[np.isfinite(pol)] = 1./2. pol[np.isfinite(pol)] = 1./2.
@@ -738,7 +738,7 @@ class overplot_radio(align_maps):
self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
self.U, self.V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.) self.U, self.V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.)
self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=1./self.vec_scale, self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=1./self.vec_scale,
scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black', label="{0:s} polarisation map".format(self.map_observer)) scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black', label="{0:s} polarization map".format(self.map_observer))
self.ax_overplot.autoscale(False) self.ax_overplot.autoscale(False)
# Display other map as contours # Display other map as contours
@@ -750,7 +750,7 @@ class overplot_radio(align_maps):
self.ax_overplot.set_xlabel(label="Right Ascension (J2000)") self.ax_overplot.set_xlabel(label="Right Ascension (J2000)")
self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1) self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1)
self.fig_overplot.suptitle("{0:s} polarisation map of {1:s} overplotted with {2:s} {3:.2f}GHz map in {4:s}.".format( self.fig_overplot.suptitle("{0:s} polarization map of {1:s} overplotted with {2:s} {3:.2f}GHz map in {4:s}.".format(
self.map_observer, obj, self.other_observer, other_freq*1e-9, self.other_unit), wrap=True) self.map_observer, obj, self.other_observer, other_freq*1e-9, self.other_unit), wrap=True)
# Display pixel scale and North direction # Display pixel scale and North direction
@@ -770,7 +770,7 @@ class overplot_radio(align_maps):
self.cr_other, = self.ax_overplot.plot(*(self.other_wcs.celestial.wcs.crpix-(1., 1.)), 'g+', transform=self.ax_overplot.get_transform(self.other_wcs)) self.cr_other, = self.ax_overplot.plot(*(self.other_wcs.celestial.wcs.crpix-(1., 1.)), 'g+', transform=self.ax_overplot.get_transform(self.other_wcs))
handles, labels = self.ax_overplot.get_legend_handles_labels() handles, labels = self.ax_overplot.get_legend_handles_labels()
handles[np.argmax([li == "{0:s} polarisation map".format(self.map_observer) for li in labels]) handles[np.argmax([li == "{0:s} polarization map".format(self.map_observer) for li in labels])
] = FancyArrowPatch((0, 0), (0, 1), arrowstyle='-', fc='w', ec='k', lw=2) ] = FancyArrowPatch((0, 0), (0, 1), arrowstyle='-', fc='w', ec='k', lw=2)
labels.append("{0:s} contour".format(self.other_observer)) labels.append("{0:s} contour".format(self.other_observer))
handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.collections[0].get_edgecolor()[0])) handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.collections[0].get_edgecolor()[0]))
@@ -829,7 +829,7 @@ class overplot_chandra(align_maps):
self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(11, 10), subplot_kw=dict(projection=self.wcs_UV)) self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(11, 10), subplot_kw=dict(projection=self.wcs_UV))
self.fig_overplot.subplots_adjust(hspace=0, wspace=0, bottom=0.1, left=0.1, top=0.8, right=1) self.fig_overplot.subplots_adjust(hspace=0, wspace=0, bottom=0.1, left=0.1, top=0.8, right=1)
# Display UV intensity map with polarisation vectors # Display UV intensity map with polarization vectors
vmin, vmax = stkI[np.isfinite(stkI)].max()/1e3*self.map_convert, stkI[np.isfinite(stkI)].max()*self.map_convert vmin, vmax = stkI[np.isfinite(stkI)].max()/1e3*self.map_convert, stkI[np.isfinite(stkI)].max()*self.map_convert
for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]:
try: try:
@@ -847,7 +847,7 @@ class overplot_chandra(align_maps):
self.cbar = self.fig_overplot.colorbar(self.im, ax=self.ax_overplot, aspect=50, shrink=0.75, pad=0.025, self.cbar = self.fig_overplot.colorbar(self.im, ax=self.ax_overplot, aspect=50, shrink=0.75, pad=0.025,
label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit)) label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit))
# Display full size polarisation vectors # Display full size polarization vectors
if vec_scale is None: if vec_scale is None:
self.vec_scale = 2. self.vec_scale = 2.
pol[np.isfinite(pol)] = 1./2. pol[np.isfinite(pol)] = 1./2.
@@ -857,7 +857,7 @@ class overplot_chandra(align_maps):
self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
self.U, self.V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.) self.U, self.V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.)
self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=1./self.vec_scale, self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=1./self.vec_scale,
scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black', label="{0:s} polarisation map".format(self.map_observer)) scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black', label="{0:s} polarization map".format(self.map_observer))
self.ax_overplot.autoscale(False) self.ax_overplot.autoscale(False)
# Display other map as contours # Display other map as contours
@@ -870,7 +870,7 @@ class overplot_chandra(align_maps):
self.ax_overplot.set_xlabel(label="Right Ascension (J2000)") self.ax_overplot.set_xlabel(label="Right Ascension (J2000)")
self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1) self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1)
self.fig_overplot.suptitle("{0:s} polarisation map of {1:s} overplotted\nwith {2:s} contour in counts.".format( self.fig_overplot.suptitle("{0:s} polarization map of {1:s} overplotted\nwith {2:s} contour in counts.".format(
self.map_observer, obj, self.other_observer), wrap=True) self.map_observer, obj, self.other_observer), wrap=True)
# Display pixel scale and North direction # Display pixel scale and North direction
@@ -889,7 +889,7 @@ class overplot_chandra(align_maps):
self.cr_map, = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix-(1., 1.)), 'r+') self.cr_map, = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix-(1., 1.)), 'r+')
self.cr_other, = self.ax_overplot.plot(*(other_wcs.celestial.wcs.crpix-(1., 1.)), 'g+', transform=self.ax_overplot.get_transform(other_wcs)) self.cr_other, = self.ax_overplot.plot(*(other_wcs.celestial.wcs.crpix-(1., 1.)), 'g+', transform=self.ax_overplot.get_transform(other_wcs))
handles, labels = self.ax_overplot.get_legend_handles_labels() handles, labels = self.ax_overplot.get_legend_handles_labels()
handles[np.argmax([li == "{0:s} polarisation map".format(self.map_observer) for li in labels]) handles[np.argmax([li == "{0:s} polarization map".format(self.map_observer) for li in labels])
] = FancyArrowPatch((0, 0), (0, 1), arrowstyle='-', fc='w', ec='k', lw=2) ] = FancyArrowPatch((0, 0), (0, 1), arrowstyle='-', fc='w', ec='k', lw=2)
labels.append("{0:s} contour in counts".format(self.other_observer)) labels.append("{0:s} contour in counts".format(self.other_observer))
handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.collections[0].get_edgecolor()[0])) handles.append(Rectangle((0, 0), 1, 1, fill=False, lw=2, ec=other_cont.collections[0].get_edgecolor()[0]))
@@ -944,7 +944,7 @@ class overplot_pol(align_maps):
self.ax_overplot.set_xlabel(label="Right Ascension (J2000)") self.ax_overplot.set_xlabel(label="Right Ascension (J2000)")
self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1) self.ax_overplot.set_ylabel(label="Declination (J2000)", labelpad=-1)
self.fig_overplot.suptitle("{0:s} observation from {1:s} overplotted with polarisation vectors and Stokes I contours from {2:s}".format( self.fig_overplot.suptitle("{0:s} observation from {1:s} overplotted with polarization vectors and Stokes I contours from {2:s}".format(
obj, self.other_observer, self.map_observer), wrap=True) obj, self.other_observer, self.map_observer), wrap=True)
# Display "other" intensity map # Display "other" intensity map
@@ -965,7 +965,7 @@ class overplot_pol(align_maps):
self.cbar = self.fig_overplot.colorbar(self.im, ax=self.ax_overplot, aspect=80, shrink=0.75, pad=0.025, self.cbar = self.fig_overplot.colorbar(self.im, ax=self.ax_overplot, aspect=80, shrink=0.75, pad=0.025,
label=r"$F_{{\lambda}}$ [{0:s}]".format(self.other_unit)) label=r"$F_{{\lambda}}$ [{0:s}]".format(self.other_unit))
# Display full size polarisation vectors # Display full size polarization vectors
if vec_scale is None: if vec_scale is None:
self.vec_scale = 2. self.vec_scale = 2.
pol[np.isfinite(pol)] = 1./2. pol[np.isfinite(pol)] = 1./2.
@@ -976,7 +976,7 @@ class overplot_pol(align_maps):
self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
self.U, self.V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.) self.U, self.V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.)
self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=px_scale/self.vec_scale, scale_units='xy', pivot='mid', self.Q = self.ax_overplot.quiver(self.X[::step_vec, ::step_vec], self.Y[::step_vec, ::step_vec], self.U[::step_vec, ::step_vec], self.V[::step_vec, ::step_vec], units='xy', angles='uv', scale=px_scale/self.vec_scale, scale_units='xy', pivot='mid',
headwidth=0., headlength=0., headaxislength=0., width=2.0, linewidth=1.0, color='white', edgecolor='black', transform=self.ax_overplot.get_transform(self.wcs_UV), label="{0:s} polarisation map".format(self.map_observer)) headwidth=0., headlength=0., headaxislength=0., width=2.0, linewidth=1.0, color='white', edgecolor='black', transform=self.ax_overplot.get_transform(self.wcs_UV), label="{0:s} polarization map".format(self.map_observer))
# Display Stokes I as contours # Display Stokes I as contours
if levels is None: if levels is None:
@@ -1009,7 +1009,7 @@ class overplot_pol(align_maps):
self.legend_title = r"{0:s} image".format(self.other_observer) self.legend_title = r"{0:s} image".format(self.other_observer)
handles, labels = self.ax_overplot.get_legend_handles_labels() handles, labels = self.ax_overplot.get_legend_handles_labels()
handles[np.argmax([li == "{0:s} polarisation map".format(self.map_observer) for li in labels]) handles[np.argmax([li == "{0:s} polarization map".format(self.map_observer) for li in labels])
] = FancyArrowPatch((0, 0), (0, 1), arrowstyle='-', fc='w', ec='k', lw=2) ] = FancyArrowPatch((0, 0), (0, 1), arrowstyle='-', fc='w', ec='k', lw=2)
labels.append("{0:s} Stokes I contour".format(self.map_observer)) labels.append("{0:s} Stokes I contour".format(self.map_observer))
handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=cont_stkI.collections[0].get_edgecolor()[0])) handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=cont_stkI.collections[0].get_edgecolor()[0]))
@@ -1358,7 +1358,7 @@ class crop_map(object):
class crop_Stokes(crop_map): class crop_Stokes(crop_map):
""" """
Class to interactively crop a polarisation map to desired Region of Interest. Class to interactively crop a polarization map to desired Region of Interest.
Inherit from crop_map. Inherit from crop_map.
""" """
@@ -1440,10 +1440,10 @@ class crop_Stokes(crop_map):
2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err) 2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err)
for dataset in self.hdul_crop: for dataset in self.hdul_crop:
dataset.header['P_int'] = (P_diluted, 'Integrated polarisation degree') dataset.header['P_int'] = (P_diluted, 'Integrated polarization degree')
dataset.header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarisation degree error') dataset.header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarization degree error')
dataset.header['PA_int'] = (PA_diluted, 'Integrated polarisation angle') dataset.header['PA_int'] = (PA_diluted, 'Integrated polarization angle')
dataset.header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarisation angle error') dataset.header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarization angle error')
self.fig.canvas.draw_idle() self.fig.canvas.draw_idle()
@property @property
@@ -1738,7 +1738,7 @@ class aperture(object):
class pol_map(object): class pol_map(object):
""" """
Class to interactively study polarisation maps. Class to interactively study polarization maps.
""" """
def __init__(self, Stokes, SNRp_cut=3., SNRi_cut=3., flux_lim=None, selection=None): def __init__(self, Stokes, SNRp_cut=3., SNRi_cut=3., flux_lim=None, selection=None):
@@ -1769,7 +1769,7 @@ class pol_map(object):
# Display selected data (Default to total flux) # Display selected data (Default to total flux)
self.display() self.display()
# Display polarisation vectors in SNR_cut # Display polarization vectors in SNR_cut
self.pol_vector() self.pol_vector()
# Display integrated values in ROI # Display integrated values in ROI
self.pol_int() self.pol_int()
@@ -2474,7 +2474,7 @@ if __name__ == "__main__":
parser = argparse.ArgumentParser(description='Interactively plot the pipeline products') parser = argparse.ArgumentParser(description='Interactively plot the pipeline products')
parser.add_argument('-f', '--file', metavar='path', required=False, help='the full or relative path to the data product', type=str, default=None) parser.add_argument('-f', '--file', metavar='path', required=False, help='the full or relative path to the data product', type=str, default=None)
parser.add_argument('-p', '--snrp', metavar='snrp_cut', required=False, help='the cut in signal-to-noise for the polarisation degree', type=float, default=3.) parser.add_argument('-p', '--snrp', metavar='snrp_cut', required=False, help='the cut in signal-to-noise for the polarization degree', type=float, default=3.)
parser.add_argument('-i', '--snri', metavar='snri_cut', required=False, help='the cut in signal-to-noise for the intensity', type=float, default=3.) parser.add_argument('-i', '--snri', metavar='snri_cut', required=False, help='the cut in signal-to-noise for the intensity', type=float, default=3.)
parser.add_argument('-l', '--lim', metavar='flux_lim', nargs=2, required=False, help='limits for the intensity map', default=None) parser.add_argument('-l', '--lim', metavar='flux_lim', nargs=2, required=False, help='limits for the intensity map', default=None)
args = parser.parse_args() args = parser.parse_args()

54
package/lib/reduction.py
View File

@@ -30,7 +30,7 @@ prototypes :
Compute Stokes parameters I, Q and U and their respective correlated errors from data_array. Compute Stokes parameters I, Q and U and their respective correlated errors from data_array.
- 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 - 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
Compute polarisation degree (in %) and angle (in degree) and their respective errors. Compute polarization degree (in %) and angle (in degree) and their respective errors.
- rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, ang, SNRi_cut) -> I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers - rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, ang, SNRi_cut) -> I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers
Rotate I, Q, U given an angle in degrees using scipy functions. Rotate I, Q, U given an angle in degrees using scipy functions.
@@ -466,7 +466,7 @@ def get_error(data_array, headers, error_array=None, data_mask=None, sub_type=No
mask = zeropad(mask_c, data[0].shape).astype(bool) mask = zeropad(mask_c, data[0].shape).astype(bool)
background = np.zeros((data.shape[0])) background = np.zeros((data.shape[0]))
# wavelength dependence of the polariser filters # wavelength dependence of the polarizer filters
# estimated to less than 1% # estimated to less than 1%
err_wav = data*0.01 err_wav = data*0.01
# difference in PSFs through each polarizers # difference in PSFs through each polarizers
@@ -963,7 +963,7 @@ def polarizer_avg(data_array, error_array, data_mask, headers, FWHM=None, scale=
pol120, err120 = smooth_data(pol120_array, err120_array, data_mask, headers120, FWHM=FWHM, scale=scale, smoothing=smoothing) pol120, err120 = smooth_data(pol120_array, err120_array, data_mask, headers120, FWHM=FWHM, scale=scale, smoothing=smoothing)
else: else:
# Sum on each polarisation filter. # Sum on each polarization filter.
pol0_t = np.sum([header['exptime'] for header in headers0]) pol0_t = np.sum([header['exptime'] for header in headers0])
pol60_t = np.sum([header['exptime'] for header in headers60]) pol60_t = np.sum([header['exptime'] for header in headers60])
pol120_t = np.sum([header['exptime'] for header in headers120]) pol120_t = np.sum([header['exptime'] for header in headers120])
@@ -1070,10 +1070,10 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale
total intensity total intensity
Q_stokes : numpy.ndarray Q_stokes : numpy.ndarray
Image (2D floats) containing the Stokes parameters accounting for Image (2D floats) containing the Stokes parameters accounting for
vertical/horizontal linear polarisation intensity vertical/horizontal linear polarization intensity
U_stokes : numpy.ndarray U_stokes : numpy.ndarray
Image (2D floats) containing the Stokes parameters accounting for Image (2D floats) containing the Stokes parameters accounting for
+45/-45deg linear polarisation intensity +45/-45deg linear polarization intensity
Stokes_cov : numpy.ndarray Stokes_cov : numpy.ndarray
Covariance matrix of the Stokes parameters I, Q, U. Covariance matrix of the Stokes parameters I, Q, U.
""" """
@@ -1236,17 +1236,17 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale
2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err) 2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err)
for header in headers: for header in headers:
header['P_int'] = (P_diluted, 'Integrated polarisation degree') header['P_int'] = (P_diluted, 'Integrated polarization degree')
header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarisation degree error') header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarization degree error')
header['PA_int'] = (PA_diluted, 'Integrated polarisation angle') header['PA_int'] = (PA_diluted, 'Integrated polarization angle')
header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarisation angle error') header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarization angle error')
return I_stokes, Q_stokes, U_stokes, Stokes_cov return I_stokes, Q_stokes, U_stokes, Stokes_cov
def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers): def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers):
""" """
Compute the polarisation degree (in %) and angle (in deg) and their Compute the polarization degree (in %) and angle (in deg) and their
respective errors from given Stokes parameters. respective errors from given Stokes parameters.
---------- ----------
Inputs: Inputs:
@@ -1255,10 +1255,10 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers):
total intensity total intensity
Q_stokes : numpy.ndarray Q_stokes : numpy.ndarray
Image (2D floats) containing the Stokes parameters accounting for Image (2D floats) containing the Stokes parameters accounting for
vertical/horizontal linear polarisation intensity vertical/horizontal linear polarization intensity
U_stokes : numpy.ndarray U_stokes : numpy.ndarray
Image (2D floats) containing the Stokes parameters accounting for Image (2D floats) containing the Stokes parameters accounting for
+45/-45deg linear polarisation intensity +45/-45deg linear polarization intensity
Stokes_cov : numpy.ndarray Stokes_cov : numpy.ndarray
Covariance matrix of the Stokes parameters I, Q, U. Covariance matrix of the Stokes parameters I, Q, U.
headers : header list headers : header list
@@ -1266,21 +1266,21 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers):
---------- ----------
Returns: Returns:
P : numpy.ndarray P : numpy.ndarray
Image (2D floats) containing the polarisation degree (in %). Image (2D floats) containing the polarization degree (in %).
debiased_P : numpy.ndarray debiased_P : numpy.ndarray
Image (2D floats) containing the debiased polarisation degree (in %). Image (2D floats) containing the debiased polarization degree (in %).
s_P : numpy.ndarray s_P : numpy.ndarray
Image (2D floats) containing the error on the polarisation degree. Image (2D floats) containing the error on the polarization degree.
s_P_P : numpy.ndarray s_P_P : numpy.ndarray
Image (2D floats) containing the Poisson noise error on the Image (2D floats) containing the Poisson noise error on the
polarisation degree. polarization degree.
PA : numpy.ndarray PA : numpy.ndarray
Image (2D floats) containing the polarisation angle. Image (2D floats) containing the polarization angle.
s_PA : numpy.ndarray s_PA : numpy.ndarray
Image (2D floats) containing the error on the polarisation angle. Image (2D floats) containing the error on the polarization angle.
s_PA_P : numpy.ndarray s_PA_P : numpy.ndarray
Image (2D floats) containing the Poisson noise error on the Image (2D floats) containing the Poisson noise error on the
polarisation angle. polarization angle.
new_headers : header list new_headers : header list
Updated list of headers corresponding to the reduced images accounting Updated list of headers corresponding to the reduced images accounting
for the new orientation angle. for the new orientation angle.
@@ -1354,10 +1354,10 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
total intensity total intensity
Q_stokes : numpy.ndarray Q_stokes : numpy.ndarray
Image (2D floats) containing the Stokes parameters accounting for Image (2D floats) containing the Stokes parameters accounting for
vertical/horizontal linear polarisation intensity vertical/horizontal linear polarization intensity
U_stokes : numpy.ndarray U_stokes : numpy.ndarray
Image (2D floats) containing the Stokes parameters accounting for Image (2D floats) containing the Stokes parameters accounting for
+45/-45deg linear polarisation intensity +45/-45deg linear polarization intensity
Stokes_cov : numpy.ndarray Stokes_cov : numpy.ndarray
Covariance matrix of the Stokes parameters I, Q, U. Covariance matrix of the Stokes parameters I, Q, U.
data_mask : numpy.ndarray data_mask : numpy.ndarray
@@ -1379,10 +1379,10 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
accounting for total intensity accounting for total intensity
new_Q_stokes : numpy.ndarray new_Q_stokes : numpy.ndarray
Rotated mage (2D floats) containing the rotated Stokes parameters Rotated mage (2D floats) containing the rotated Stokes parameters
accounting for vertical/horizontal linear polarisation intensity accounting for vertical/horizontal linear polarization intensity
new_U_stokes : numpy.ndarray new_U_stokes : numpy.ndarray
Rotated image (2D floats) containing the rotated Stokes parameters Rotated image (2D floats) containing the rotated Stokes parameters
accounting for +45/-45deg linear polarisation intensity. accounting for +45/-45deg linear polarization intensity.
new_Stokes_cov : numpy.ndarray new_Stokes_cov : numpy.ndarray
Updated covariance matrix of the Stokes parameters I, Q, U. Updated covariance matrix of the Stokes parameters I, Q, U.
new_headers : header list new_headers : header list
@@ -1497,10 +1497,10 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err) 2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err)
for header in new_headers: for header in new_headers:
header['P_int'] = (P_diluted, 'Integrated polarisation degree') header['P_int'] = (P_diluted, 'Integrated polarization degree')
header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarisation degree error') header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarization degree error')
header['PA_int'] = (PA_diluted, 'Integrated polarisation angle') header['PA_int'] = (PA_diluted, 'Integrated polarization angle')
header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarisation angle error') header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarization angle error')
return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_data_mask, new_headers return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_data_mask, new_headers

18
package/overplot_MRK463E.py
View File

@@ -10,17 +10,11 @@ Stokes_Xr = fits.open("./data/MRK463E/Chandra/X_ray_crop.fits")
levels = np.geomspace(1., 99., 7) levels = np.geomspace(1., 99., 7)
# A = overplot_chandra(Stokes_UV, Stokes_Xr) A = overplot_chandra(Stokes_UV, Stokes_Xr, norm=LogNorm())
# A.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=20.0, zoom=1, savename='./plots/MRK463E/Chandra_overplot.pdf') A.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=5, zoom=1, savename='./plots/MRK463E/Chandra_overplot.pdf')
A.write_to(path1="./data/MRK463E/FOC_data_Chandra.fits", path2="./data/MRK463E/Chandra_data.fits", suffix="aligned")
B = overplot_chandra(Stokes_UV, Stokes_Xr, norm=LogNorm())
B.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=5, zoom=1, savename='./plots/MRK463E/Chandra_overplot_forced.pdf')
B.write_to(path1="./data/MRK463E/FOC_data_Chandra.fits", path2="./data/MRK463E/Chandra_data.fits", suffix="aligned")
# C = overplot_pol(Stokes_UV, Stokes_IR)
# C.plot(SNRp_cut=3.0, SNRi_cut=20.0, savename='./plots/MRK463E/IR_overplot.pdf')
levels = np.array([0.8, 2, 5, 10, 20, 50])/100.*Stokes_UV[0].header['photflam'] levels = np.array([0.8, 2, 5, 10, 20, 50])/100.*Stokes_UV[0].header['photflam']
D = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm()) B = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm())
D.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=5, norm=LogNorm(8.5e-18, 2.5e-15), savename='./plots/MRK463E/IR_overplot_forced.pdf') B.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=5, norm=LogNorm(8.5e-18, 2.5e-15), savename='./plots/MRK463E/IR_overplot.pdf')
D.write_to(path1="./data/MRK463E/FOC_data_WFPC.fits", path2="./data/MRK463E/WFPC_data.fits", suffix="aligned") B.write_to(path1="./data/MRK463E/FOC_data_WFPC.fits", path2="./data/MRK463E/WFPC_data.fits", suffix="aligned")