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

add 180° as a convention to PA

Browse Source
This commit is contained in:
Tibeuleu
2022-11-09 16:05:17 +01:00
parent 359e5b401f
commit 47f0ead97c
32 changed files with 36 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/FOC_reduction.py
View File

@@ -128,7 +128,7 @@ def main():
# Data binning
rebin = True
if rebin:
pxsize = 0.05
pxsize = 0.10
px_scale = 'arcsec' #pixel, arcsec or full
rebin_operation = 'sum' #sum or average
# Alignement
@@ -136,7 +136,7 @@ def main():
display_data = False
# Smoothing
smoothing_function = 'combine' #gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine
smoothing_FWHM = 0.10 #If None, no smoothing is done
smoothing_FWHM = 0.20 #If None, no smoothing is done
smoothing_scale = 'arcsec' #pixel or arcsec
# Rotation
rotate_stokes = True #rotation to North convention can give erroneous results
@@ -146,7 +146,7 @@ def main():
final_display = True
# Polarization map output
figname = 'NGC1068_FOC' #target/intrument name
figtype = '_combine_FWHM010' #additionnal informations
figtype = '_combine_FWHM020' #additionnal informations
SNRp_cut = 5. #P measurments with SNR>3
SNRi_cut = 50. #I measurments with SNR>30, which implies an uncertainty in P of 4.7%.
step_vec = 1 #plot all vectors in the array. if step_vec = 2, then every other vector will be plotted

14
src/comparison_Kishimoto.py
View File

@@ -1,5 +1,5 @@
#!/usr/bin/env python
from lib.reduction import align_data, princ_angle
from lib.reduction import align_data, crop_array, princ_angle
from lib.deconvolve import zeropad
from matplotlib.colors import LogNorm
from os.path import join as path_join
@@ -46,7 +46,7 @@ for d in [data_S, data_K]:
d[i][np.isnan(d[i])] = 0.
d['P'] = np.where(d['I']>0.,np.sqrt(d['Q']**2+d['U']**2)/d['I'],0.)
d['sP'] = np.where(d['I']>0.,np.sqrt((d['Q']**2*d['sQ']**2+d['U']**2*d['sU']**2)/(d['Q']**2+d['U']**2)+((d['Q']/d['I'])**2+(d['U']/d['I'])**2)*d['sI']**2)/d['I'],0.)
d['PA'] = (90./np.pi)*np.arctan2(d['U'],d['Q'])
d['PA'] = princ_angle((90./np.pi)*np.arctan2(d['U'],d['Q'])+180.)
d['SNRp'] = np.zeros(d['P'].shape)
d['SNRp'][d['sP']>0.] = d['P'][d['sP']>0.]/d['sP'][d['sP']>0.]
d['SNRi'] = np.zeros(d['I'].shape)
@@ -65,7 +65,11 @@ ax = fig.add_subplot(111, projection=wcs)
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.12, 0.01, 0.75])
im0 = ax.imshow(data_S['I']*convert_flux,norm=LogNorm(data_S['I'][data_S['I']>0].min()*convert_flux,data_S['I'][data_S['I']>0].max()*convert_flux),origin='lower',cmap='gray',label=r"$I_{STOKES}$ through this pipeline")
#im0 = ax.imshow(data_S['I']*convert_flux,norm=LogNorm(data_S['I'][data_S['I']>0].min()*convert_flux,data_S['I'][data_S['I']>0].max()*convert_flux),origin='lower',cmap='gray',label=r"$I_{STOKES}$ through this pipeline")
#im0 = ax.imshow(data_K['P']*100.,vmin=0.,vmax=100.,origin='lower',cmap='inferno',label=r"$P$ through Kishimoto's pipeline")
im0 = ax.imshow(data_S['P']*100.,vmin=0.,vmax=100.,origin='lower',cmap='inferno',label=r"$P$ through this pipeline")
#im0 = ax.imshow(data_K['PA'],vmin=0.,vmax=360.,origin='lower',cmap='inferno',label=r"$\theta_P$ through Kishimoto's pipeline")
#im0 = ax.imshow(data_S['PA'],vmin=0.,vmax=360.,origin='lower',cmap='inferno',label=r"$\theta_P$ through this pipeline")
quiv0 = ax.quiver(data_S['X'],data_S['Y'],data_S['xy_U'],data_S['xy_V'],units='xy',angles='uv',scale=0.5,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,color='b',alpha=0.75, label="PA through this pipeline")
quiv1 = ax.quiver(data_K['X'],data_K['Y'],data_K['xy_U'],data_K['xy_V'],units='xy',angles='uv',scale=0.5,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,color='r',alpha=0.75, label="PA through Kishimoto's pipeline")
@@ -79,7 +83,9 @@ ax.coords[1].set_axislabel_position('l')
ax.coords[1].set_ticklabel_position('l')
#ax.axis('equal')
cbar = plt.colorbar(im0, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
#cbar = plt.colorbar(im0, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
cbar = plt.colorbar(im0, cax=cbar_ax, label=r"$P$ [%]")
#cbar = plt.colorbar(im0, cax=cbar_ax, label=r"$\theta_P$ [°]")
plt.rcParams.update({'font.size': 15})
ax.legend(loc='upper right')

26
src/lib/plots.py
View File

@@ -450,20 +450,16 @@ class align_maps(object):
self.other_map = other_map
self.wcs_map = deepcopy(WCS(self.map[0])).celestial
# if self.wcs_map.naxis == 4:
# self.wcs_map = WCS(self.map[0],naxis=[1,2]).deepcopy()
# self.map[0].data = self.map[0].data[0,0]
# elif self.wcs_map.naxis == 3:
# self.wcs_map = WCS(self.map[0],naxis=[1,2]).deepcopy()
# self.map[0].data = self.map[0].data[1]
if len(self.map[0].data.shape) == 4:
self.map[0].data = self.map[0].data[0,0]
elif len(self.map[0].data.shape) == 3:
self.map[0].data = self.map[0].data[1]
self.wcs_other = deepcopy(WCS(self.other_map[0])).celestial
# if self.wcs_other.naxis == 4:
# self.wcs_other = WCS(self.other_map[0],naxis=[1,2]).deepcopy()
# self.other_map[0].data = self.other_map[0].data[0,0]
# elif self.wcs_other.naxis == 3:
# self.wcs_other = WCS(self.other_map[0],naxis=[1,2]).deepcopy()
# self.other_map[0].data = self.other_map[0].data[1]
if len(self.other_map[0].data.shape) == 4:
self.other_map[0].data = self.other_map[0].data[0,0]
elif len(self.other_map[0].data.shape) == 3:
self.other_map[0].data = self.other_map[0].data[1]
try:
convert_flux = self.map[0].header['photflam']
@@ -1837,7 +1833,7 @@ class pol_map(object):
P_cut = np.sqrt(Q_cut**2+U_cut**2)/I_cut
P_cut_err = np.sqrt((Q_cut**2*Q_cut_err**2 + U_cut**2*U_cut_err**2 + 2.*Q_cut*U_cut*QU_cut_err)/(Q_cut**2 + U_cut**2) + ((Q_cut/I_cut)**2 + (U_cut/I_cut)**2)*I_cut_err**2 - 2.*(Q_cut/I_cut)*IQ_cut_err - 2.*(U_cut/I_cut)*IU_cut_err)/I_cut
PA_cut = princ_angle(np.degrees((1./2.)*np.arctan2(U_cut,Q_cut)))
PA_cut = princ_angle(np.degrees((1./2.)*np.arctan2(U_cut,Q_cut))+180.)
PA_cut_err = princ_angle(np.degrees((1./(2.*(Q_cut**2+U_cut**2)))*np.sqrt(U_cut**2*Q_cut_err**2 + Q_cut**2*U_cut_err**2 - 2.*Q_cut*U_cut*QU_cut_err)))
else:
@@ -1862,7 +1858,7 @@ class pol_map(object):
P_reg = np.sqrt(Q_reg**2+U_reg**2)/I_reg
P_reg_err = np.sqrt((Q_reg**2*Q_reg_err**2 + U_reg**2*U_reg_err**2 + 2.*Q_reg*U_reg*QU_reg_err)/(Q_reg**2 + U_reg**2) + ((Q_reg/I_reg)**2 + (U_reg/I_reg)**2)*I_reg_err**2 - 2.*(Q_reg/I_reg)*IQ_reg_err - 2.*(U_reg/I_reg)*IU_reg_err)/I_reg
PA_reg = princ_angle((90./np.pi)*np.arctan2(U_reg,Q_reg))
PA_reg = princ_angle((90./np.pi)*np.arctan2(U_reg,Q_reg)+180.)
PA_reg_err = (90./(np.pi*(Q_reg**2+U_reg**2)))*np.sqrt(U_reg**2*Q_reg_err**2 + Q_reg**2*U_reg_err**2 - 2.*Q_reg*U_reg*QU_reg_err)
new_cut = np.logical_and(self.region, self.cut)
@@ -1879,7 +1875,7 @@ class pol_map(object):
P_cut = np.sqrt(Q_cut**2+U_cut**2)/I_cut
P_cut_err = np.sqrt((Q_cut**2*Q_cut_err**2 + U_cut**2*U_cut_err**2 + 2.*Q_cut*U_cut*QU_cut_err)/(Q_cut**2 + U_cut**2) + ((Q_cut/I_cut)**2 + (U_cut/I_cut)**2)*I_cut_err**2 - 2.*(Q_cut/I_cut)*IQ_cut_err - 2.*(U_cut/I_cut)*IU_cut_err)/I_cut
PA_cut = princ_angle((90./np.pi)*np.arctan2(U_cut,Q_cut))
PA_cut = princ_angle((90./np.pi)*np.arctan2(U_cut,Q_cut)+180.)
PA_cut_err = (90./(np.pi*(Q_cut**2+U_cut**2)))*np.sqrt(U_cut**2*Q_cut_err**2 + Q_cut**2*U_cut_err**2 - 2.*Q_cut*U_cut*QU_cut_err)
if hasattr(self, 'cont'):

6
src/lib/reduction.py
View File

@@ -1310,7 +1310,7 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
P_diluted = np.sqrt(Q_diluted**2+U_diluted**2)/I_diluted
P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted**2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err)
PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted,Q_diluted))
PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted,Q_diluted)+180.)
PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err**2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err)
for header in headers:
@@ -1370,7 +1370,7 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers):
P = np.zeros(I_stokes.shape)
P[mask] = I_pol[mask]/I_stokes[mask]
PA = np.zeros(I_stokes.shape)
PA[mask] = princ_angle((90./np.pi)*np.arctan2(U_stokes[mask],Q_stokes[mask]))
PA[mask] = princ_angle((90./np.pi)*np.arctan2(U_stokes[mask],Q_stokes[mask])+180.)
if (P>1).any():
print("WARNING : found {0:d} pixels for which P > 1".format(P[P>1.].size))
@@ -1569,7 +1569,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
P_diluted = np.sqrt(Q_diluted**2+U_diluted**2)/I_diluted
P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted**2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err)
PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted,Q_diluted))
PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted,Q_diluted)+180.)
PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err**2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err)
for header in new_headers:

18
src/overplot.py
View File

@@ -6,35 +6,35 @@ from copy import deepcopy
from lib.plots import overplot_radio, overplot_pol, align_pol
from matplotlib.colors import LogNorm
#Stokes_UV = fits.open("../data/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.fits")
Stokes_UV = fits.open("../data/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.fits")
#Stokes_18GHz = fits.open("../data/IC5063_x3nl030/radio/IC5063.18GHz.fits")
#Stokes_24GHz = fits.open("../data/IC5063_x3nl030/radio/IC5063.24GHz.fits")
#Stokes_103GHz = fits.open("../data/IC5063_x3nl030/radio/I5063_103GHz.fits")
#Stokes_229GHz = fits.open("../data/IC5063_x3nl030/radio/I5063_229GHz.fits")
#Stokes_357GHz = fits.open("../data/IC5063_x3nl030/radio/I5063_357GHz.fits")
#Stokes_S2 = fits.open("../data/IC5063_x3nl030/POLARIZATION_COMPARISON/S2_rot_crop.fits")
#Stokes_IR = fits.open("../data/IC5063_x3nl030/IR/u2e65g01t_c0f_rot.fits")
Stokes_IR = fits.open("../data/IC5063_x3nl030/IR/u2e65g01t_c0f_rot.fits")
#levelsMorganti = np.array([1.,2.,3.,8.,16.,32.,64.,128.])
#
##levels18GHz = np.array([0.6, 1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_18GHz[0].data.max()
#levels18GHz = levelsMorganti*0.28*1e-3
#A = overplot_radio(Stokes_UV, Stokes_18GHz)
#A.plot(levels=levels18GHz, SNRp_cut=3.0, SNRi_cut=60.0, savename='../plots/IC5063_x3nl030/18GHz_overplot_forced.png')
#
##levels24GHz = np.array([1.,1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_24GHz[0].data.max()
#levels24GHz = levelsMorganti*0.46*1e-3
#B = overplot_radio(Stokes_UV, Stokes_24GHz)
#B.plot(levels=levels24GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/24GHz_overplot_forced.png')
#
#levels103GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_103GHz[0].data[Stokes_103GHz[0].data > 0.]))
#C = overplot_radio(Stokes_UV, Stokes_103GHz)
#C.plot(levels=levels103GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/103GHz_overplot_forced.png')
#
#levels229GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_229GHz[0].data[Stokes_229GHz[0].data > 0.]))
#D = overplot_radio(Stokes_UV, Stokes_229GHz)
#D.plot(levels=levels229GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/229GHz_overplot_forced.png')
#
#levels357GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_357GHz[0].data[Stokes_357GHz[0].data > 0.]))
#E = overplot_radio(Stokes_UV, Stokes_357GHz)
#E.plot(levels=levels357GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/357GHz_overplot_forced.png')
@@ -42,8 +42,8 @@ from matplotlib.colors import LogNorm
#F = overplot_pol(Stokes_UV, Stokes_S2)
#F.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/S2_overplot_forced.png', norm=LogNorm(vmin=5e-20,vmax=5e-18))
#G = overplot_pol(Stokes_UV, Stokes_IR, cmap='inferno')
#G.plot(SNRp_cut=3.0, SNRi_cut=60.0, savename='../plots/IC5063_x3nl030/IR_overplot_forced.png', norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r')
G = overplot_pol(Stokes_UV, Stokes_IR, cmap='inferno')
G.plot(SNRp_cut=3.0, SNRi_cut=60.0, savename='../plots/IC5063_x3nl030/IR_overplot_forced.png', norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r')
#data_folder1 = "../data/M87/POS1/"
#plots_folder1 = "../plots/M87/POS1/"