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

new images for corrected polarization angle

Browse Source
This commit is contained in:
Tibeuleu
2022-12-15 14:43:19 +01:00
parent 93f43394e2
commit 6684b3ab49
120 changed files with 37 additions and 31 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
src/FOC_reduction.py
View File

@@ -124,11 +124,11 @@ def main():
display_crop = False
# Error estimation
error_sub_shape = (15,15)
display_error = True
display_error = False
# 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,17 +136,17 @@ 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
rotate_data = False #rotation to North convention can give erroneous results
# Final crop
crop = False #Crop to desired ROI
final_display = True
final_display = False
# 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

4
src/comparison_Kishimoto.py
View File

@@ -20,7 +20,7 @@ data_K = {}
data_S = {}
for d,i in zip(['I','Q','U','P','PA','sI','sQ','sU','sP','sPA'],[0,1,2,5,8,(3,0,0),(3,1,1),(3,2,2),6,9]):
data_K[d] = np.loadtxt(path_join(root_dir_K,d+'.txt'))
with fits.open(path_join(root_dir_data_S,'NGC1068_old_FOC_bin10px.fits')) as f:
with fits.open(path_join(root_dir_data_S,'NGC1068_FOC_bin10px.fits')) as f:
if not type(i) is int:
data_S[d] = np.sqrt(f[i[0]].data[i[1],i[2]])
else:
@@ -110,7 +110,7 @@ print("This pipeline : average sI/I={0:.2f} ; sQ/Q={1:.2f} ; sU/U={2:.2f} ; sP/P
print("Kishimoto's pipeline : average sI/I={0:.2f} ; sQ/Q={1:.2f} ; sU/U={2:.2f} ; sP/P={3:.2f}".format(np.mean(data_K['sI'][data_S['mask']]/data_K['I'][data_S['mask']]),np.mean(data_K['sQ'][data_S['mask']]/data_K['Q'][data_S['mask']]),np.mean(data_K['sU'][data_S['mask']]/data_K['U'][data_S['mask']]),np.mean(data_K['sP'][data_S['mask']]/data_K['P'][data_S['mask']])))
for d,i in zip(['I','Q','U','P','PA','sI','sQ','sU','sP','sPA'],[0,1,2,5,8,(3,0,0),(3,1,1),(3,2,2),6,9]):
data_K[d] = np.loadtxt(path_join(root_dir_K,d+'.txt'))
with fits.open(path_join(root_dir_data_S,'NGC1068_K_FOC_bin10px.fits')) as f:
with fits.open(path_join(root_dir_data_S,'NGC1068_FOC_bin10px.fits')) as f:
if not type(i) is int:
data_S[d] = np.sqrt(f[i[0]].data[i[1],i[2]])
else:

24
src/lib/plots.py
View File

@@ -59,8 +59,8 @@ def princ_angle(ang):
A = np.array(ang)
while np.any(A < 0.):
A[A<0.] = A[A<0.]+360.
while np.any(A >= 360.):
A[A>=360.] = A[A>=360.]-360.
while np.any(A >= 180.):
A[A>=180.] = A[A>=180.]-180.
if type(ang) == type(A):
return A
else:
@@ -336,7 +336,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
elif display.lower() in ['pa','pang','pol_ang']:
# Display polarization degree map
display='pa'
vmin, vmax = 0., 360.
vmin, vmax = 0., 180.
im = ax.imshow(princ_angle(pang.data), vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.)
cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\theta_P$ [°]")
elif display.lower() in ['s_p','pol_err','pol_deg_err']:
@@ -674,6 +674,9 @@ class overplot_radio(align_maps):
north_dir = AnchoredDirectionArrows(self.ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color='w', arrow_props={'ec': None, 'fc': 'w', 'alpha': 1,'lw': 2})
self.ax.add_artist(north_dir)
self.cr_map, = self.ax.plot(*self.wcs_map.wcs.crpix, 'r+')
crpix_other = self.wcs_map.world_to_pixel(self.wcs_other.pixel_to_world(*self.wcs_other.wcs.crpix))
self.cr_other, = self.ax.plot(*crpix_other, 'g+')
if not(savename is None):
self.fig2.savefig(savename,bbox_inches='tight',dpi=200)
@@ -728,7 +731,7 @@ class overplot_pol(align_maps):
#Display Stokes I as contours
levels_stkI = np.rint(np.linspace(10,99,10))/100.*np.max(stkI.data[stkI.data > 0.]*convert_flux)
cont_stkI = self.ax.contour(stkI.data*convert_flux, transform=self.ax.get_transform(self.wcs_UV), levels=levels_stkI, colors='grey', alpha=0.)
cont_stkI = self.ax.contour(stkI.data*convert_flux, transform=self.ax.get_transform(self.wcs_UV), levels=levels_stkI, colors='grey', alpha=0.5)
#self.ax.clabel(cont_stkI, inline=True, fontsize=8)
self.ax.autoscale(False)
@@ -752,6 +755,8 @@ class overplot_pol(align_maps):
cbar_ax = self.fig2.add_axes([0.95, 0.12, 0.01, 0.75])
cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
self.ax.set(xlabel="Right Ascension (J2000)", ylabel="Declination (J2000)", title="{0:s} overplotted with polarization vectors and Stokes I contours from HST/FOC".format(obj))
#Display pixel scale and North direction
fontprops = fm.FontProperties(size=16)
px_size = self.wcs_UV.wcs.get_cdelt()[0]*3600.
@@ -760,8 +765,9 @@ class overplot_pol(align_maps):
north_dir = AnchoredDirectionArrows(self.ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color='w', arrow_props={'ec': None, 'fc': 'w', 'alpha': 1,'lw': 2})
self.ax.add_artist(north_dir)
self.ax.set(xlabel="Right Ascension (J2000)", ylabel="Declination (J2000)", title="{0:s} overplotted with polarization vectors and Stokes I contours from HST/FOC".format(obj))
self.cr_map, = self.ax.plot(*self.wcs_map.wcs.crpix, 'r+')
crpix_other = self.wcs_map.world_to_pixel(self.wcs_other.pixel_to_world(*self.wcs_other.wcs.crpix))
self.cr_other, = self.ax.plot(*crpix_other, 'g+')
if not(savename is None):
self.fig2.savefig(savename,bbox_inches='tight',dpi=200)
@@ -1744,7 +1750,7 @@ class pol_map(object):
label = r"$P$ [%]"
elif self.display_selection.lower() in ['pol_ang']:
self.data = princ_angle(self.PA)
vmin, vmax = 0, 360.
vmin, vmax = 0, 180.
label = r"$\theta_{P}$ [°]"
elif self.display_selection.lower() in ['snri']:
s_I = np.sqrt(self.IQU_cov[0,0])
@@ -1833,7 +1839,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 = 360.-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)))
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:
@@ -1858,7 +1864,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 = 360.-princ_angle((90./np.pi)*np.arctan2(U_reg,Q_reg))
PA_reg = princ_angle((90./np.pi)*np.arctan2(U_reg,Q_reg))
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)

14
src/lib/reduction.py
View File

@@ -80,8 +80,8 @@ def princ_angle(ang):
A = np.array(ang)
while np.any(A < 0.):
A[A<0.] = A[A<0.]+360.
while np.any(A >= 360.):
A[A>=360.] = A[A>=360.]-360.
while np.any(A >= 180.):
A[A>=180.] = A[A>=180.]-180.
if type(ang) == type(A):
return A
else:
@@ -480,7 +480,6 @@ def get_error_hist(data_array, headers, error_array=None, data_mask=None,
date_time = np.array([datetime.strptime(d,'%Y-%m-%d;%H:%M:%S')
for d in date_time])
for i, image in enumerate(data):
filt_obs[headers[i]['filtnam1']] += 1
#Compute the Count-rate histogram for the image
n_mask = np.logical_and(mask,image>0.)
@@ -499,9 +498,10 @@ def get_error_hist(data_array, headers, error_array=None, data_mask=None,
#bkg = np.percentile(image[image<hist_max],25.)
#bkg = 0.95*hist_max
if display:
filt_obs[headers[i]['filtnam1']] += 1
ax_h.plot(bin_centers,hist,'+',color="C{0:d}".format(i),alpha=0.8,label=headers[i]['filtnam1']+' (Obs '+str(filt_obs[headers[i]['filtnam1']])+')')
ax_h.plot([bkg,bkg],[hist.min(), hist.max()],'x--',color="C{0:d}".format(i),alpha=0.8)
print(headers[i]['filtnam1']+' ('+str(date_time[i])+') : n_bins =',n_bins,'; bkg = {0:.2e}'.format(bkg))
#print(headers[i]['filtnam1']+' ('+str(date_time[i])+') : n_bins =',n_bins,'; bkg = {0:.2e}'.format(bkg))
error_bkg[i] *= bkg
# Quadratically add uncertainties in the "correction factors" (see Kishimoto 1999)
@@ -993,7 +993,7 @@ def align_data(data_array, headers, error_array=None, background=None,
raise ValueError("All images in data_array must have same shape as\
ref_data")
if (error_array is None) or (background is None):
_, error_array, headers, background = get_error(data_array, headers, return_background=True)
_, error_array, headers, background = get_error_hist(data_array, headers, return_background=True)
# Crop out any null edges
#(ref_data must be cropped as well)
@@ -1517,7 +1517,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 = 360.-princ_angle((90./np.pi)*np.arctan2(U_diluted,Q_diluted))
PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted,Q_diluted))
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:
@@ -1776,7 +1776,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 = 360.-princ_angle((90./np.pi)*np.arctan2(U_diluted,Q_diluted))
PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted,Q_diluted))
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:

16
src/overplot.py
View File

@@ -6,7 +6,7 @@ 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_old_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")
@@ -20,30 +20,30 @@ 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=30.0, savename='../plots/IC5063_x3nl030/18GHz_overplot_old.png')
A.plot(levels=levels18GHz, SNRp_cut=2.0, SNRi_cut=15.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=30.0, savename='../plots/IC5063_x3nl030/24GHz_overplot_old.png')
B.plot(levels=levels24GHz, SNRp_cut=2.0, SNRi_cut=15.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=30.0, savename='../plots/IC5063_x3nl030/103GHz_overplot_old.png')
C.plot(levels=levels103GHz, SNRp_cut=2.0, SNRi_cut=15.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=30.0, savename='../plots/IC5063_x3nl030/229GHz_overplot_old.png')
D.plot(levels=levels229GHz, SNRp_cut=2.0, SNRi_cut=15.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=30.0, savename='../plots/IC5063_x3nl030/357GHz_overplot_old.png')
E.plot(levels=levels357GHz, SNRp_cut=2.0, SNRi_cut=15.0, savename='../plots/IC5063_x3nl030/357GHz_overplot_forced.png')
#F = overplot_pol(Stokes_UV, Stokes_S2)
#F.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/S2_overplot_old.png', norm=LogNorm(vmin=5e-20,vmax=5e-18))
#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=30.0, savename='../plots/IC5063_x3nl030/IR_overplot_old.png', norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r')
G.plot(SNRp_cut=2.0, SNRi_cut=15.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/"