diff --git a/plots/IC5063_x3nl030/IC5063_FOC_center_image.png b/plots/IC5063_x3nl030/IC5063_FOC_center_image.png index 0041d38..f37dc85 100755 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_center_image.png and b/plots/IC5063_x3nl030/IC5063_FOC_center_image.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png new file mode 100644 index 0000000..67c0959 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png new file mode 100644 index 0000000..ce3e1ab Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png new file mode 100644 index 0000000..a4a28af Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png new file mode 100644 index 0000000..c895051 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png new file mode 100644 index 0000000..95d156b Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png new file mode 100644 index 0000000..5142d9b Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png new file mode 100644 index 0000000..d31e6a1 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png new file mode 100644 index 0000000..f652092 Binary files /dev/null and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020.png new file mode 100644 index 0000000..126bb77 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_IQU.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_IQU.png new file mode 100644 index 0000000..ce5c55a Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_IQU.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_I_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_I_err.png new file mode 100644 index 0000000..ae19cc4 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_I_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P.png new file mode 100644 index 0000000..09dddfe Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_err.png new file mode 100644 index 0000000..24d32f7 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_flux.png new file mode 100644 index 0000000..1c5fc27 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_P_flux.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRi.png new file mode 100644 index 0000000..c9e3ec7 Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRp.png new file mode 100644 index 0000000..f3c9fcd Binary files /dev/null and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_SNRp.png differ diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 31b5b22..d04b73a 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -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', @@ -81,13 +81,13 @@ def main(): # infiles = ['x25d0401t_c0f.fits','x25d0402t_c0f.fits','x25d0403t_c0f.fits'] # globals()['plots_folder'] = "../plots/PictorA_x25d040/" - globals()['data_folder'] = "../data/3C273_x0u20/" - infiles = ['x0u20101t_c0f.fits','x0u20102t_c0f.fits','x0u20103t_c0f.fits','x0u20104t_c0f.fits','x0u20105t_c0f.fits','x0u20106t_c0f.fits','x0u20201t_c0f.fits','x0u20202t_c0f.fits','x0u20203t_c0f.fits','x0u20204t_c0f.fits','x0u20205t_c0f.fits','x0u20206t_c0f.fits','x0u20301t_c0f.fits','x0u20302t_c0f.fits','x0u20303t_c0f.fits','x0u20304t_c0f.fits','x0u20305t_c0f.fits','x0u20306t_c0f.fits'] - globals()['plots_folder'] = "../plots/3C273_x0u20/" +# globals()['data_folder'] = "../data/3C273_x0u20/" +# infiles = ['x0u20101t_c0f.fits','x0u20102t_c0f.fits','x0u20103t_c0f.fits','x0u20104t_c0f.fits','x0u20105t_c0f.fits','x0u20106t_c0f.fits','x0u20201t_c0f.fits','x0u20202t_c0f.fits','x0u20203t_c0f.fits','x0u20204t_c0f.fits','x0u20205t_c0f.fits','x0u20206t_c0f.fits','x0u20301t_c0f.fits','x0u20302t_c0f.fits','x0u20303t_c0f.fits','x0u20304t_c0f.fits','x0u20305t_c0f.fits','x0u20306t_c0f.fits'] +# globals()['plots_folder'] = "../plots/3C273_x0u20/" ## Reduction parameters # Deconvolution - deconvolve = True + deconvolve = False if deconvolve: psf = 'gaussian' #Can be user-defined as well psf_FWHM = 0.15 @@ -97,12 +97,12 @@ def main(): # Cropping display_crop = False # Error estimation - error_sub_shape = (75,75) + error_sub_shape = (150,150) display_error = False # Data binning rebin = True if rebin: - pxsize = 0.15 + pxsize = 0.05 px_scale = 'arcsec' #pixel or arcsec rebin_operation = 'sum' #sum or average # Alignement @@ -110,16 +110,16 @@ def main(): display_data = False # Smoothing smoothing_function = 'combine' #gaussian_after, gaussian or combine - smoothing_FWHM = 0.15 #If None, no smoothing is done + smoothing_FWHM = 0.10 #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 # Polarization map output - figname = '3C273_FOC' #target/intrument name - figtype = '_combine_FWHM015_deconvolved' #additionnal informations - SNRp_cut = 5. #P measurments with SNR>3 - SNRi_cut = 20. #I measurments with SNR>30, which implies an uncertainty in P of 4.7%. + figname = 'NGC1068_FOC' #target/intrument name + figtype = '_combine_FWHM020' #additionnal informations + SNRp_cut = 50. #P measurments with SNR>3 + SNRi_cut = 350. #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 diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 3021a28..8469063 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -626,7 +626,7 @@ def rebin_array(data_array, error_array, headers, pxsize, scale, def align_data(data_array, headers, error_array=None, upsample_factor=1., - ref_data=None, ref_center=None, return_shifts=True): + ref_data=None, ref_center=None, return_shifts=False): """ Align images in data_array using cross correlation, and rescale them to wider images able to contain any rotation of the reference image. @@ -851,7 +851,7 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1., g_rc = np.array([np.exp(-0.5*(dist_rc/stdev)**2),]*data_array.shape[0]) # Apply weighted combination smoothed[r,c] = (1.-data_mask[r,c])*np.sum(data_array*weight*g_rc)/np.sum(weight*g_rc) - error[r,c] = np.sqrt(np.sum(weight*g_rc**2))/np.sum(weight*g_rc) + error[r,c] = (1.-data_mask[r,c])*np.sqrt(np.sum(weight*g_rc**2))/np.sum(weight*g_rc) # Nan handling error[np.isnan(smoothed)] = 0. @@ -871,7 +871,7 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1., with warnings.catch_warnings(record=True) as w: g_rc = np.exp(-0.5*(dist_rc/stdev)**2)/(2.*np.pi*stdev**2) smoothed[i][r,c] = (1.-data_mask[r,c])*np.sum(image*g_rc) - error[i][r,c] = np.sqrt(np.sum(error_array[i]*g_rc**2)) + error[i][r,c] = (1.-data_mask[r,c])*np.sqrt(np.sum(error_array[i]*g_rc**2)) # Nan handling error[i][np.isnan(smoothed[i])] = 0. @@ -1440,3 +1440,72 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, new_Stokes_cov[np.isnan(new_Stokes_cov)] = fmax return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_data_mask, new_headers + +def rotate_data(data_array, error_array, data_mask, headers, ang): + """ + Use scipy.ndimage.rotate to rotate I_stokes to an angle, and a rotation + matrix to rotate Q, U of a given angle in degrees and update header + orientation keyword. + ---------- + Inputs: + data_array : numpy.ndarray + Array of images (2D floats) to be rotated by angle ang. + error_array : numpy.ndarray + Array of error associated to images in data_array. + headers : header list + List of headers corresponding to the reduced images. + ang : float + Rotation angle (in degrees) that should be applied to the Stokes + parameters + ---------- + Returns: + new_data_array : numpy.ndarray + Updated array containing the rotated images. + new_error_array : numpy.ndarray + Updated array containing the rotated errors. + new_headers : header list + Updated list of headers corresponding to the reduced images accounting + for the new orientation angle. + """ + #Rotate I_stokes, Q_stokes, U_stokes using rotation matrix + alpha = ang*np.pi/180. + + #Rotate original images using scipy.ndimage.rotate + new_data_array = [] + new_error_array = [] + for i in range(data_array.shape[0]): + new_data_array.append(sc_rotate(data_array[i], ang, reshape=False, + cval=0.)) + new_error_array.append(sc_rotate(error_array[i], ang, reshape=False, + cval=error_array.mean())) + new_data_array = np.array(new_data_array) + new_data_mask = sc_rotate(data_mask, ang, reshape=False, cval=True) + new_error_array = np.array(new_error_array) + + for i in range(new_data_array.shape[0]): + new_data_array[i][new_data_array[i] < 0.] = 0. + + #Update headers to new angle + new_headers = [] + mrot = np.array([[np.cos(-alpha), -np.sin(-alpha)], + [np.sin(-alpha), np.cos(-alpha)]]) + for header in headers: + new_header = deepcopy(header) + new_header['orientat'] = header['orientat'] + ang + + new_wcs = WCS(header).deepcopy() + if new_wcs.wcs.has_cd(): # CD matrix + del new_wcs.wcs.cd + keys = ['CD1_1','CD1_2','CD2_1','CD2_2'] + for key in keys: + new_header.remove(key, ignore_missing=True) + new_wcs.wcs.cdelt = 3600.*np.sqrt(np.sum(new_wcs.wcs.get_pc()**2,axis=1)) + elif new_wcs.wcs.has_pc(): # PC matrix + CDELT + newpc = np.dot(mrot, new_wcs.wcs.get_pc()) + new_wcs.wcs.pc = newpc + new_wcs.wcs.set() + new_header.update(new_wcs.to_header()) + + new_headers.append(new_header) + + return new_data_array, new_error_array, new_data_mask, new_headers