reduce NGC1068 for comparison with Kishimoto's pipeline

src/FOC_reduction.py

@@ -18,12 +18,13 @@ from astropy.wcs import WCS
 
 ##### User inputs
 ## Input and output locations
-#globals()['data_folder'] = "../data/NGC1068_x274020/"
-#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']
@@ -90,13 +91,13 @@ from astropy.wcs import WCS
 #globals()['plots_folder'] = "../plots/3C273_x0u20/"
 
 #BEWARE: 5 observations separated by 1 year each (1995, 1996, 1997, 1998, 1999)
-globals()['data_folder'] = "../data/M87/POS1/"
+#globals()['data_folder'] = "../data/M87/POS1/"
 #globals()['infiles'] = ['x2py010ct_c0f.fits','x2py010dt_c0f.fits','x2py010et_c0f.fits','x2py010ft_c0f.fits'] #1995
 #globals()['infiles'] = ['x3be010ct_c0f.fits','x3be010dt_c0f.fits','x3be010et_c0f.fits','x3be010ft_c0f.fits'] #1996
 #globals()['infiles'] = ['x43r010km_c0f.fits','x43r010mm_c0f.fits','x43r010om_c0f.fits','x43r010rm_c0f.fits'] #1997
 #globals()['infiles'] = ['x43r110kr_c0f.fits','x43r110mr_c0f.fits','x43r110or_c0f.fits','x43r110rr_c0f.fits'] #1998
-globals()['infiles'] = ['x43r210kr_c0f.fits','x43r210mr_c0f.fits','x43r210or_c0f.fits','x43r210rr_c0f.fits'] #1999
-globals()['plots_folder'] = "../plots/M87/POS1/"
+#globals()['infiles'] = ['x43r210kr_c0f.fits','x43r210mr_c0f.fits','x43r210or_c0f.fits','x43r210rr_c0f.fits'] #1999
+#globals()['plots_folder'] = "../plots/M87/POS1/"
 
 #BEWARE: 5 observations separated by 1 year each (1995, 1996, 1997, 1998, 1999)
 #globals()['data_folder'] = "../data/M87/POS3/"
@@ -128,28 +129,28 @@ def main():
     # Data binning
     rebin = True
     if rebin:
-        pxsize = 0.05
-        px_scale = 'arcsec'         #pixel, arcsec or full
+        pxsize = 10
+        px_scale = 'pixel'         #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 = 0.10           #If None, no smoothing is done
+    smoothing_FWHM = None           #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 = 'M87_POS1_1999_FOC'         #target/intrument name
-    figtype = '_combine_FWHM010'    #additionnal informations
+    figname = 'NGC1068_K_FOC'         #target/intrument name
+    figtype = '_bin10px'    #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 = 0    #plot all vectors in the array. if step_vec = 2, then every other vector will be plotted
+    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
 
     ##### Pipeline start

src/lib/plots.py

@@ -1703,12 +1703,14 @@ class pol_map(object):
         ax.add_artist(self.north_dir)
 
     def display(self, fig=None, ax=None):
+        norm = None
         if self.display_selection is None:
             self.display_selection = "total_flux"
 
         if self.display_selection.lower() in ['total_flux']:
             self.data = self.I*self.convert_flux
-            vmin, vmax = 0., np.max(self.data[self.data > 0.])
+            vmin, vmax = np.min(self.data[self.cut])/10., np.max(self.data[self.data > 0.])
+            norm = LogNorm(vmin, vmax)
             label = r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]"
         elif self.display_selection.lower() in ['pol_flux']:
             self.data = self.I*self.convert_flux*self.P
@@ -1738,10 +1740,16 @@ class pol_map(object):
                 ax = self.ax
             if hasattr(self, 'im'):
                 self.im.remove()
+            if not norm is None:
+                self.im = ax.imshow(self.data, norm=norm, aspect='equal', cmap='inferno')
+            else:
                 self.im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno')
             self.cbar = plt.colorbar(self.im, cax=self.cbar_ax, label=label)
             fig.canvas.draw_idle()
             return self.im
+        else:
+            if not norm is None:
+                im = ax.imshow(self.data, norm=norm, aspect='equal', cmap='inferno')
             else:
                 im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno')
             ax.set_xlim(0,self.data.shape[1])

src/lib/reduction.py

@@ -1475,24 +1475,23 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
     new_U_stokes = np.zeros(shape)
     new_Stokes_cov = np.zeros((*Stokes_cov.shape[:-2],*shape))
 
-    for i in range(shape[0]):
-        for j in range(shape[1]):
-            new_I_stokes[i,j], new_Q_stokes[i,j], new_U_stokes[i,j] = np.dot(mrot, np.array([I_stokes[i,j], Q_stokes[i,j], U_stokes[i,j]])).T
-            new_Stokes_cov[:,:,i,j] = np.dot(mrot, np.dot(Stokes_cov[:,:,i,j], mrot.T))
-
     #Rotate original images using scipy.ndimage.rotate
-    new_I_stokes = sc_rotate(new_I_stokes, ang, order=1, reshape=False, cval=0.)
-    new_Q_stokes = sc_rotate(new_Q_stokes, ang, order=1, reshape=False, cval=0.)
-    new_U_stokes = sc_rotate(new_U_stokes, ang, order=1, reshape=False, cval=0.)
+    new_I_stokes = sc_rotate(I_stokes, ang, order=1, reshape=False, cval=0.)
+    new_Q_stokes = sc_rotate(Q_stokes, ang, order=1, reshape=False, cval=0.)
+    new_U_stokes = sc_rotate(U_stokes, ang, order=1, reshape=False, cval=0.)
     new_data_mask = sc_rotate(data_mask.astype(float)*10., ang, order=1, reshape=False, cval=0.)
     new_data_mask[new_data_mask < 2] = 0.
     new_data_mask = new_data_mask.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, order=1,
-                    reshape=False, cval=0.)
+            new_Stokes_cov[i,j] = sc_rotate(Stokes_cov[i,j], ang, order=1, reshape=False, cval=0.)
         new_Stokes_cov[i,i] = np.abs(new_Stokes_cov[i,i])
 
+    for i in range(shape[0]):
+        for j in range(shape[1]):
+            new_I_stokes[i,j], new_Q_stokes[i,j], new_U_stokes[i,j] = np.dot(mrot, np.array([new_I_stokes[i,j], new_Q_stokes[i,j], new_U_stokes[i,j]])).T
+            new_Stokes_cov[:,:,i,j] = np.dot(mrot, np.dot(new_Stokes_cov[:,:,i,j], mrot.T))
+
     #Update headers to new angle
     new_headers = []
     mrot = np.array([[np.cos(-alpha), -np.sin(-alpha)],