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

redo images without axis error (still buggy)

Browse Source
This commit is contained in:
Thibault Barnouin
2021-07-01 20:44:38 +02:00
parent 7958246655
commit 4e0861129d
113 changed files with 40 additions and 30 deletions
Download Patch File Download Diff File Expand all files Collapse all files

32
src/FOC_reduction.py
View File

@@ -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',
@@ -77,6 +77,10 @@ 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/"
## Reduction parameters
# Deconvolution
deconvolve = False
@@ -108,10 +112,10 @@ 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 = 'NGC1068_FOC' #target/intrument name
figname = 'PictorA_FOC' #target/intrument name
figtype = '_combine_FWHM020_rot' #additionnal informations
SNRp_cut = 3. #P measurments with SNR>3
SNRi_cut = 30 #I measurments with SNR>30, which implies an uncertainty in P of 4.7%.
SNRp_cut = 20. #P measurments with SNR>3
SNRi_cut = 120. #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
##### Pipeline start
@@ -193,12 +197,12 @@ def main():
## Step 5:
# Plot polarization map (Background is either total Flux, Polarization degree or Polarization degree error).
proj_plots.polarization_map(copy.deepcopy(Stokes_test), rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype, plots_folder=plots_folder, display=None)
proj_plots.polarization_map(copy.deepcopy(Stokes_test), rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P_flux", plots_folder=plots_folder, display='Pol_Flux')
proj_plots.polarization_map(copy.deepcopy(Stokes_test), rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P", plots_folder=plots_folder, display='Pol_deg')
proj_plots.polarization_map(copy.deepcopy(Stokes_test), rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P_err", plots_folder=plots_folder, display='Pol_deg_err')
proj_plots.polarization_map(copy.deepcopy(Stokes_test), rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_SNRi", plots_folder=plots_folder, display='SNRi')
proj_plots.polarization_map(copy.deepcopy(Stokes_test), rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_SNRp", plots_folder=plots_folder, display='SNRp')
proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype, plots_folder=plots_folder, display=None)
proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P_flux", plots_folder=plots_folder, display='Pol_Flux')
proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P", plots_folder=plots_folder, display='Pol_deg')
proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_P_err", plots_folder=plots_folder, display='Pol_deg_err')
proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_SNRi", plots_folder=plots_folder, display='SNRi')
proj_plots.polarization_map(copy.deepcopy(Stokes_test), data_mask, rectangle=None, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, step_vec=step_vec, savename=figname+figtype+"_SNRp", plots_folder=plots_folder, display='SNRp')
return 0

25
src/lib/plots.py
View File

@@ -167,7 +167,7 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""):
return 0
def polarization_map(Stokes, rectangle=None, SNRp_cut=3., SNRi_cut=30., step_vec=1,
def polarization_map(Stokes, data_mask, rectangle=None, SNRp_cut=3., SNRi_cut=30., step_vec=1,
savename=None, plots_folder="", display=None):
"""
Plots polarization map from Stokes HDUList.
@@ -235,6 +235,7 @@ def polarization_map(Stokes, rectangle=None, SNRp_cut=3., SNRi_cut=30., step_vec
SNRi[np.isnan(SNRi)] = 0.
pol.data[SNRi < SNRi_cut] = np.nan
data_mask = (1.-data_mask).astype(bool)
mask = (SNRp > SNRp_cut) * (SNRi > SNRi_cut)
# Look for pixel of max polarization
@@ -340,16 +341,16 @@ def polarization_map(Stokes, rectangle=None, SNRp_cut=3., SNRi_cut=30., step_vec
PA_int_err = (90./(np.pi*(Q_int**2 + U_int**2)))*np.sqrt(U_int**2*Q_int_err**2 + Q_int**2*U_int_err**2 - 2.*Q_int*U_int*QU_int_err)
# Compute integrated parameters and associated errors for all pixels
n_pix = stkI.data.size
I_diluted = stkI.data.sum()
Q_diluted = stkQ.data.sum()
U_diluted = stkU.data.sum()
I_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[0,0]))
Q_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[1,1]))
U_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[2,2]))
IQ_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[0,1]**2))
IU_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[0,2]**2))
QU_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[1,2]**2))
n_pix = stkI.data[data_mask].size
I_diluted = stkI.data[data_mask].sum()
Q_diluted = stkQ.data[data_mask].sum()
U_diluted = stkU.data[data_mask].sum()
I_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[0,0][data_mask]))
Q_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[1,1][data_mask]))
U_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[2,2][data_mask]))
IQ_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[0,1][data_mask]**2))
IU_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[0,2][data_mask]**2))
QU_diluted_err = np.sqrt(n_pix)*np.sqrt(np.sum(stk_cov.data[1,2][data_mask]**2))
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)
@@ -359,7 +360,7 @@ def polarization_map(Stokes, rectangle=None, SNRp_cut=3., SNRi_cut=30., step_vec
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)
#PA_diluted_err = P_diluted_err/(2.*P_diluted)*180./np.pi
ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(pivot_wav,sci_not(I_int*convert_flux,I_int_err*convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_int*100.,P_int_err*100.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_int,PA_int_err), color='white', fontsize=16, xy=(0.01, 0.92), xycoords='axes fraction')
ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(pivot_wav,sci_not(I_diluted*convert_flux,I_diluted_err*convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_diluted*100.,P_diluted_err*100.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_diluted,PA_diluted_err), color='white', fontsize=16, xy=(0.01, 0.92), xycoords='axes fraction')
ax.coords.grid(True, color='white', ls='dotted', alpha=0.5)
ax.coords[0].set_axislabel('Right Ascension (J2000)')

8
src/lib/reduction.py
View File

@@ -1153,9 +1153,9 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
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 + A*coeff_stokes[2,1]*(np.sin(2.*theta[1]*Q_stokes) - np.cos(2.*theta[1]*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 + A*coeff_stokes[2,2]*(np.sin(2.*theta[2]*Q_stokes) - np.cos(2.*theta[2]*U_stokes)))
Stokes_cov[0,0] += (dI_dtheta1**2*sigma_theta[0]**2 + dI_dtheta2**2*sigma_theta[1]**2 + dI_dtheta3**2*sigma_theta[2]**2)
Stokes_cov[1,1] += (dQ_dtheta1**2*sigma_theta[0]**2 + dQ_dtheta2**2*sigma_theta[1]**2 + dQ_dtheta3**2*sigma_theta[2]**2)
Stokes_cov[2,2] += (dU_dtheta1**2*sigma_theta[0]**2 + dU_dtheta2**2*sigma_theta[1]**2 + dU_dtheta3**2*sigma_theta[2]**2)
#Stokes_cov[0,0] += (dI_dtheta1**2*sigma_theta[0]**2 + dI_dtheta2**2*sigma_theta[1]**2 + dI_dtheta3**2*sigma_theta[2]**2)
#Stokes_cov[1,1] += (dQ_dtheta1**2*sigma_theta[0]**2 + dQ_dtheta2**2*sigma_theta[1]**2 + dQ_dtheta3**2*sigma_theta[2]**2)
#Stokes_cov[2,2] += (dU_dtheta1**2*sigma_theta[0]**2 + dU_dtheta2**2*sigma_theta[1]**2 + dU_dtheta3**2*sigma_theta[2]**2)
plt.imshow(np.abs(Stokes_cov[0,0]/I_stokes)*100., origin='lower')
plt.colorbar()
plt.show()
@@ -1337,7 +1337,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
new_I_stokes = sc_rotate(new_I_stokes, ang, reshape=False, cval=0.)
new_Q_stokes = sc_rotate(new_Q_stokes, ang, reshape=False, cval=0.)
new_U_stokes = sc_rotate(new_U_stokes, ang, reshape=False, cval=0.)
new_data_mask = sc_rotate(data_mask, ang, reshape=False, cval=True)
new_data_mask = (sc_rotate(data_mask.astype(float)*10., ang, reshape=False, cval=10.).astype(int)).astype(bool)
for i in range(3):
for j in range(3):
new_Stokes_cov[i,j] = sc_rotate(new_Stokes_cov[i,j], ang,