correction on principal angle computation

src/FOC_reduction.py

@@ -18,18 +18,18 @@ from astropy.wcs import WCS
 
 ##### User inputs
 ## Input and output locations
-globals()['data_folder'] = "../data/NGC1068_x274020/"
-#globals()['infiles'] = ['xn1c400.fits','xn2c400.fits','xn3c400.fits']
-globals()['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']
-#psf_file = 'NGC1068_f253m00.fits'
-globals()['plots_folder'] = "../plots/NGC1068_x274020/"
+#globals()['data_folder'] = "../data/NGC1068_x274020/"
+##globals()['infiles'] = ['xn1c400.fits','xn2c400.fits','xn3c400.fits']
+#globals()['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']
+##psf_file = 'NGC1068_f253m00.fits'
+#globals()['plots_folder'] = "../plots/NGC1068_x274020/"
 
-#globals()['data_folder'] = "../data/IC5063_x3nl030/"
-#globals()['infiles'] = ['x3nl0301r_c0f.fits','x3nl0302r_c0f.fits','x3nl0303r_c0f.fits']
-##psf_file = 'IC5063_f502m00.fits'
-#globals()['plots_folder'] = "../plots/IC5063_x3nl030/"
+globals()['data_folder'] = "../data/IC5063_x3nl030/"
+globals()['infiles'] = ['x3nl0301r_c0f.fits','x3nl0302r_c0f.fits','x3nl0303r_c0f.fits']
+#psf_file = 'IC5063_f502m00.fits'
+globals()['plots_folder'] = "../plots/IC5063_x3nl030/"
 
 #globals()['data_folder'] = "../data/NGC1068_x14w010/"
 #globals()['infiles'] = ['x14w0101t_c0f.fits','x14w0102t_c0f.fits','x14w0103t_c0f.fits',
@@ -129,26 +129,26 @@ def main():
     # Data binning
     rebin = True
     if rebin:
-        pxsize = 10
-        px_scale = 'pixel'         #pixel, arcsec or full
+        pxsize = 0.10
+        px_scale = 'arcsec'         #pixel, arcsec or full
         rebin_operation = 'sum'     #sum or average
     # Alignement
     align_center = 'image'          #If None will align image to image center
     display_data = False
     # Smoothing
     smoothing_function = 'combine'  #gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine
-    smoothing_FWHM = None           #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 = False
+    final_display = True
     # Polarization map output
-    figname = 'NGC1068_K_FOC'         #target/intrument name
-    figtype = '_bin10px'    #additionnal informations
-    SNRp_cut = 5.    #P measurments with SNR>3
+    figname = 'IC5063_FOC'         #target/intrument name
+    figtype = '_combine_FWHM020'    #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%.
     step_vec = 1    #plot all vectors in the array. if step_vec = 2, then every other vector will be plotted
                     # if step_vec = 0 then all vectors are displayed at full length

src/comparison_Kishimoto.py

@@ -80,12 +80,8 @@ print('From my pipeline :\n', "P = {0:.2f} ± {1:.2f} %\n".format(data_S['P_dil'
 print("From Kishimoto's pipeline :\n", "P = {0:.2f} ± {1:.2f} %\n".format(data_K['P_dil']*100.,data_K['sP_dil']*100.), "PA = {0:.2f} ± {1:.2f} °".format(data_K['PA_dil'],data_K['sPA_dil']))
 
 #compare different types of error
-xx, yy = np.indices(data_S['mask'].shape)
-mask_ind = np.array([[y,x] for y,x in zip(yy[data_S['mask']],xx[data_S['mask']])])
-index = mask_ind[np.random.randint(len(mask_ind))]
-print("My pipeline : sI/I={0:.2f} ; sQ/Q={1:.2f} ; sU/U={2:.2f} ; sP/P={3:.2f}".format(np.mean(data_S['sI'][index[0],index[1]]/data_S['I'][index[0],index[1]]),np.mean(data_S['sQ'][index[0],index[1]]/data_S['Q'][index[0],index[1]]),np.mean(data_S['sU'][index[0],index[1]]/data_S['U'][index[0],index[1]]),np.mean(data_S['sP'][index[0],index[1]]/data_S['P'][index[0],index[1]])))
-print("Kishimoto's pipeline : sI/I={0:.2f} ; sQ/Q={1:.2f} ; sU/U={2:.2f} ; sP/P={3:.2f}".format(np.mean(data_K['sI'][index[0],index[1]]/data_K['I'][index[0],index[1]]),np.mean(data_K['sQ'][index[0],index[1]]/data_K['Q'][index[0],index[1]]),np.mean(data_K['sU'][index[0],index[1]]/data_K['U'][index[0],index[1]]),np.mean(data_K['sP'][index[0],index[1]]/data_K['P'][index[0],index[1]])))
-print("For random pixel in cut at {}".format(index))
+print("My pipeline : average sI/I={0:.2f} ; sQ/Q={1:.2f} ; sU/U={2:.2f} ; sP/P={3:.2f}".format(np.mean(data_S['sI'][data_S['mask']]/data_S['I'][data_S['mask']]),np.mean(data_S['sQ'][data_S['mask']]/data_S['Q'][data_S['mask']]),np.mean(data_S['sU'][data_S['mask']]/data_S['U'][data_S['mask']]),np.mean(data_S['sP'][data_S['mask']]/data_S['P'][data_S['mask']])))
+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:

src/lib/fits.py

@@ -13,6 +13,7 @@ import numpy as np
 from astropy.io import fits
 from astropy import wcs
 from lib.convex_hull import image_hull, clean_ROI
+from lib.plots import princ_angle
 import matplotlib.pyplot as plt
 
 
@@ -66,6 +67,7 @@ def get_obs_data(infiles, data_folder="", compute_flux=False):
             new_wcs.wcs.cdelt = new_cdelt
             for key, val in new_wcs.to_header().items():
                 header[key] = val
+        header['orientat'] = princ_angle(float(header['orientat']))
 
     if compute_flux:
         for i in range(len(infiles)):

src/lib/plots.py

@@ -51,12 +51,12 @@ from astropy.io import fits
 
 def princ_angle(ang):
     """
-    Return the principal angle in the 0-180° quadrant.
+    Return the principal angle in the -180° to 180° quadrant.
     """
-    while ang < 0.:
-        ang += 180.
+    while ang <= -180.:
+        ang += 360.
     while ang > 180.:
-        ang -= 180.
+        ang -= 360.
     return ang
 
 
@@ -1808,8 +1808,8 @@ 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_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)
+            PA_reg = np.degrees((1./2.)*np.arctan2(U_reg,Q_reg))
+            PA_reg_err = princ_angle(np.degrees((1./(2.*(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)))
 
         if hasattr(self, 'cont'):
             for coll in self.cont.collections:

src/lib/reduction.py

@@ -488,11 +488,13 @@ def get_error(data_array, headers, error_array=None, data_mask=None,
         rectangle.append([minima[1], minima[0], sub_shape[1], sub_shape[0], 0., 'r'])
         # Compute error : root mean square of the background
         sub_image = image[minima[0]:minima[0]+sub_shape[0],minima[1]:minima[1]+sub_shape[1]]
-        #error =  np.std(sub_image)    # Previously computed using standard deviation over the background
+        #bkg =  np.std(sub_image)    # Previously computed using standard deviation over the background
         bkg = np.sqrt(np.sum((sub_image-sub_image.mean())**2)/sub_image.size)
         error_bkg[i] *= bkg
-
+        
+        #Substract background
         data_array[i] = np.abs(data_array[i] - sub_image.mean())
+        
         # Quadratically add uncertainties in the "correction factors" (see Kishimoto 1999)
         #wavelength dependence of the polariser filters
         #estimated to less than 1%
@@ -1087,7 +1089,7 @@ def polarizer_avg(data_array, error_array, data_mask, headers, FWHM=None,
             if not(FWHM is None) and (smoothing.lower() in ['gaussian','gauss']):
                 # Smooth by convoluting with a gaussian each polX image.
                 pol_array, polerr_array = smooth_data(pol_array, polerr_array,
-                        data_mask, headers_array, FWHM=FWHM, scale=scale)
+                        data_mask, headers, FWHM=FWHM, scale=scale)
                 pol0, pol60, pol120 = pol_array
                 err0, err60, err120 = polerr_array
 
@@ -1294,8 +1296,8 @@ 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_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 = np.degrees((1./2.)*np.arctan2(U_diluted,Q_diluted))
+        PA_diluted_err = princ_angle(np.degrees((1./(2.*(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:
             header['P_int'] = (P_diluted, 'Integrated polarization degree')
@@ -1470,7 +1472,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
         for i,head in enumerate(headers):
             ang[i] = -head['orientat']
         ang = ang.mean()
-    alpha = ang*np.pi/180.
+    alpha = np.radians(ang)
     mrot = np.array([[1., 0., 0.],
                     [0., np.cos(2.*alpha), np.sin(2.*alpha)],
                     [0, -np.sin(2.*alpha), np.cos(2.*alpha)]])
@@ -1553,8 +1555,8 @@ 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_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 = np.degrees((1./2.)*np.arctan2(U_diluted,Q_diluted))
+    PA_diluted_err = princ_angle(np.degrees((1./(1.*(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:
         header['P_int'] = (P_diluted, 'Integrated polarization degree')