From 3cc9769eb127deaa751f79ec8cd8f82d47ac9122 Mon Sep 17 00:00:00 2001 From: enloro Date: Thu, 3 Jun 2021 09:03:56 -0700 Subject: [PATCH] rotation step and change smoothin at the end of reduction --- src/FOC_reduction_ELR.py | 72 +++++++++++++++++++++------------------- 1 file changed, 37 insertions(+), 35 deletions(-) diff --git a/src/FOC_reduction_ELR.py b/src/FOC_reduction_ELR.py index 2037c38..9563d08 100644 --- a/src/FOC_reduction_ELR.py +++ b/src/FOC_reduction_ELR.py @@ -98,10 +98,11 @@ infiles = {'data_folder':'../data/NGC1068_x274020/',\ 'x274020ht.c0f.fits','x274020it.c0f.fits']} config = {'conserve_flux':True,'upsample_factor':1,\ - 'downsample_factor':8} + 'downsample_factor':8,\ + 'rotate':True,'reshape':True} fig = {'cmap':'magma','vmin':-20,'vmax':200,'font_size':40.0,'label_size':20,'lw':3,\ - 'SNRp_cut':3,'SNRi_cut':4,'scalevec':0.05,'step_vec':1,\ + 'SNRp_cut':1,'SNRi_cut':4,'scalevec':0.05,'step_vec':1,\ 'vec_legend':10,'figname':'NGC1068_FOC_ELR'} ##################### @@ -116,42 +117,47 @@ data_array_c = cross_corr(data_array) #resample image print('--- Downsample image') data_array_cr_m, data_array_cr_s = rebin_array(data_array_c,config) -#smoothing -print('--- Smoothing') -kernel = Gaussian2DKernel(x_stddev=1) -data_array_crs_m = [] -data_array_crs_s = [] -for ii in range(len(data_array_cr_m)): - data_array_crs_m.append(convolve(data_array_cr_m[ii], kernel)) - data_array_crs_s.append(convolve(data_array_cr_m[ii], kernel)) - +#compute stokes +print('--- Combine HWP PAs') #combine HWP PAs -pol0 = (data_array_crs_m[0]+data_array_crs_m[1]+data_array_crs_m[2]) -pol60 = (data_array_crs_m[3]+data_array_crs_m[4]+data_array_crs_m[5]+data_array_crs_m[6]) -pol120 = (data_array_crs_m[7]+data_array_crs_m[8]) - -#Stokes parameters +pol0 = (data_array_cr_m[0]+data_array_cr_m[1]+data_array_cr_m[2]) +pol60 = (data_array_cr_m[3]+data_array_cr_m[4]+data_array_cr_m[5]+data_array_cr_m[6]) +pol120 = (data_array_cr_m[7]+data_array_cr_m[8]) +#Stokes +print('--- Compute Stokes') I_stokes = (2./3.)*(pol0 + pol60 + pol120) Q_stokes = (2./3.)*(2*pol0 - pol60 - pol120) U_stokes = (2./np.sqrt(3.))*(pol60 - pol120) - -#rotate Stokes -#PA = 46.81 #deg -#I_stokes = ndimage.rotate(I_stokes,-PA,reshape=True) -#Q_stokes = ndimage.rotate(Q_stokes,-PA,reshape=True) -#U_stokes = ndimage.rotate(U_stokes,-PA,reshape=True) -#Q_stokes = Q_stokes*np.cos(np.deg2rad(PA)) - U_stokes*np.sin(np.deg2rad(PA)) -#U_stokes = Q_stokes*np.sin(np.deg2rad(PA)) + U_stokes*np.cos(np.deg2rad(PA)) - +#Rotate Stokes +if config['rotate']: + print('--- Rotate Stokes') + ima = fits.open(infiles['data_folder']+infiles['filenames'][0]) + PA = ima[0].header['ORIENTAT'] + #roate images + I_stokes = ndimage.rotate(I_stokes,-PA,reshape=config['reshape']) + Q_stokes = ndimage.rotate(Q_stokes,-PA,reshape=config['reshape']) + U_stokes = ndimage.rotate(U_stokes,-PA,reshape=config['reshape']) + #rotate stokes + Q_stokes = Q_stokes*np.cos(np.deg2rad(PA)) - U_stokes*np.sin(np.deg2rad(PA)) + U_stokes = Q_stokes*np.sin(np.deg2rad(PA)) + U_stokes*np.cos(np.deg2rad(PA)) +#smoothing +print('--- Smoothing') +kernel = Gaussian2DKernel(x_stddev=1) +I_stokes = convolve(I_stokes,kernel) +Q_stokes = convolve(Q_stokes,kernel) +U_stokes = convolve(U_stokes,kernel) +#remove annoying pixels due to rotation +I_stokes[I_stokes == 0] = np.nan +Q_stokes[Q_stokes == 0] = np.nan +U_stokes[U_stokes == 0] = np.nan #Polarimetry PI = np.sqrt(Q_stokes*Q_stokes + U_stokes*U_stokes) P = PI/I_stokes*100 PA = 0.5*arctan2(U_stokes,Q_stokes)*180./np.pi+90 -s_P = np.sqrt(2.)*(I_stokes)**(-0.5) +s_P = (np.sqrt(2.)*(I_stokes)**(-0.5)) ### Correct this line!!! s_PA = s_P/(P/100.)*180./np.pi - ### STEP 6. image to FITS with updated WCS img = fits.open(infiles['data_folder']+infiles['filenames'][0]) new_wcs = WCS(naxis=2) @@ -162,7 +168,6 @@ new_wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"] new_wcs.wcs.cdelt = [img[0].header['CD1_1']*config['downsample_factor'],\ img[0].header['CD1_2']*config['downsample_factor']] -#hdu_ori = WCS(img[0]) stkI=fits.PrimaryHDU(data=I_stokes,header=new_wcs.to_header()) pol=fits.PrimaryHDU(data=P,header=new_wcs.to_header()) pang=fits.PrimaryHDU(data=PA,header=new_wcs.to_header()) @@ -175,19 +180,16 @@ pang_err=fits.PrimaryHDU(data=s_PA,header=new_wcs.to_header()) pxscale = stkI.header['CDELT1'] #apply quality cuts -SNRp = pol.data/pol_err.data -pol.data[SNRp < fig['SNRp_cut']] = np.nan +pol.data[stkI.data < 10] = np.nan -SNRi = stkI.data/np.std(stkI.data[0:10,0:10]) -SNRi = fits.PrimaryHDU(SNRi,header=stkI.header) -pol.data[SNRi.data < fig['SNRi_cut']] = np.nan -levelsI = 2**(np.arange(2,20,0.5)) +#levelsI = 2**(np.arange(2,20,0.5)) +levelsI = 2**(np.arange(3,20,0.5)) fig1 = plt.figure(figsize=(11,10)) gc = FITSFigure(stkI,figure=fig1) gc.show_colorscale(cmap=fig['cmap'],vmin=fig['vmin'],vmax=fig['vmax']) -gc.show_contour(SNRi,levels=levelsI,colors='grey',linewidths=0.5) +gc.show_contour(levels=levelsI,colors='grey',linewidths=0.5) gc.show_vectors(pol,pang,scale=fig['scalevec'],\ step=fig['step_vec'],color='black',linewidth=2.0) gc.show_vectors(pol,pang,scale=fig['scalevec'],\