remove axis error computation on polarization components

src/FOC_reduction.py

@@ -105,7 +105,7 @@ def main():
     align_center = 'image'        #If None will align image to image center
     display_data = False
     # Smoothing
-    smoothing_function = 'combine'  #gaussian_after, gaussian or combine
+    smoothing_function = 'gaussian'  #gaussian_after, gaussian or combine
     smoothing_FWHM = 0.20           #If None, no smoothing is done
     smoothing_scale = 'arcsec'       #pixel or arcsec
     # Rotation
@@ -113,7 +113,7 @@ def main():
     rotate_data = False              #rotation to North convention can give erroneous results
     # Polarization map output
     figname = 'NGC1068_FOC'         #target/intrument name
-    figtype = '_combine_FWHM020_waeP'    #additionnal informations
+    figtype = '_combine_FWHM020_wae'    #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
@@ -182,7 +182,7 @@ def main():
     # FWHM of FOC have been estimated at about 0.03" across 1500-5000 Angstrom band, which is about 2 detector pixels wide
     # see Jedrzejewski, R.; Nota, A.; Hack, W. J., A Comparison Between FOC and WFPC2
     # Bibcode : 1995chst.conf...10J
-    I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta = proj_red.compute_Stokes(data_array, error_array, data_mask, headers, FWHM=smoothing_FWHM, scale=smoothing_scale, smoothing=smoothing_function)
+    I_stokes, Q_stokes, U_stokes, Stokes_cov = proj_red.compute_Stokes(data_array, error_array, data_mask, headers, FWHM=smoothing_FWHM, scale=smoothing_scale, smoothing=smoothing_function)
 
     ## Step 3:
     # Rotate images to have North up
@@ -193,9 +193,9 @@ def main():
             [np.sin(-alpha), np.cos(-alpha)]])
         rectangle[0:2] = np.dot(mrot, np.asarray(rectangle[0:2]))+np.array(data_array.shape[1:])/2
         rectangle[4] = alpha
-        I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta, headers, data_mask = proj_red.rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta, data_mask, headers, -ref_header['orientat'], SNRi_cut=None)
+        I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, data_mask = proj_red.rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, -ref_header['orientat'], SNRi_cut=None)
     # Compute polarimetric parameters (polarization degree and angle).
-    P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta, headers)
+    P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers)
 
     ## Step 4:
     # crop to desired region of interest (roi)

src/lib/reduction.py

@@ -66,16 +66,6 @@ globals()['theta'] = np.array([180.*np.pi/180., 60.*np.pi/180., 120.*np.pi/180.]
 # Uncertainties on the orientation of the polarizers' axes taken to be 3deg (see Nota et. al 1996, p36; Robinson & Thomson 1995)
 globals()['sigma_theta'] = np.array([3.*np.pi/180., 3.*np.pi/180., 3.*np.pi/180.])
 
-def princ_angle(ang):
-    """
-    Return the principal angle in the 0-180° quadrant.
-    """
-    while ang < 0.:
-        ang += 180.
-    while ang > 180.:
-        ang -= 180.
-    return ang
-
 
 def get_row_compressor(old_dimension, new_dimension, operation='sum'):
     """
@@ -1068,7 +1058,7 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
     # Routine for the FOC instrument
     if instr == 'FOC':
         # Get image from each polarizer and covariance matrix
-        pol_array, pol_cov_m = polarizer_avg(data_array, error_array, data_mask,
+        pol_array, pol_cov = polarizer_avg(data_array, error_array, data_mask,
                 headers, FWHM=FWHM, scale=scale, smoothing=smoothing)
         pol0, pol60, pol120 = pol_array
 
@@ -1103,15 +1093,14 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
 
         # Orientation and error for each polarizer
         fmax = np.finfo(np.float64).max
-        pol_flux = 2.*np.array([pol0/transmit[0], pol60/transmit[1], pol120/transmit[2]])#np.array([pol0, pol60, pol120])#
-        pol_cov = np.array([[4.*pol_cov_m[i,j]/(transmit[i]*transmit[j]) for j in range(pol_cov_m.shape[1])] for i in range(pol_cov_m.shape[0])])
+        pol_flux = np.array([pol0, pol60, pol120])
 
         coeff_stokes = np.zeros((3,3))
         # Coefficients linking each polarizer flux to each Stokes parameter
         for i in range(3):
-            coeff_stokes[0,i] = pol_eff[(i+1)%3]*pol_eff[(i+2)%3]*np.sin(-2.*theta[(i+1)%3]+2.*theta[(i+2)%3])#*2./transmit[i]
-            coeff_stokes[1,i] = (-pol_eff[(i+1)%3]*np.sin(2.*theta[(i+1)%3]) + pol_eff[(i+2)%3]*np.sin(2.*theta[(i+2)%3]))#*2./transmit[i]
-            coeff_stokes[2,i] = (pol_eff[(i+1)%3]*np.cos(2.*theta[(i+1)%3]) - pol_eff[(i+2)%3]*np.cos(2.*theta[(i+2)%3]))#*2./transmit[i]
+            coeff_stokes[0,i] = pol_eff[(i+1)%3]*pol_eff[(i+2)%3]*np.sin(-2.*theta[(i+1)%3]+2.*theta[(i+2)%3])*2./transmit[i]
+            coeff_stokes[1,i] = (-pol_eff[(i+1)%3]*np.sin(2.*theta[(i+1)%3]) + pol_eff[(i+2)%3]*np.sin(2.*theta[(i+2)%3]))*2./transmit[i]
+            coeff_stokes[2,i] = (pol_eff[(i+1)%3]*np.cos(2.*theta[(i+1)%3]) - pol_eff[(i+2)%3]*np.cos(2.*theta[(i+2)%3]))*2./transmit[i]
 
         # Normalization parameter for Stokes parameters computation
         A = coeff_stokes[0,:].sum()
@@ -1157,24 +1146,6 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
         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)
         dU_dtheta = np.array([dU_dtheta1, dU_dtheta2, dU_dtheta3])
         
-        # Compute the derivative of each polarization component with respect to the polarizer orientation
-        mask = I_stokes > 0.
-        dP_dtheta1 = np.zeros(I_stokes.shape)*fmax
-        dP_dtheta2 = np.zeros(I_stokes.shape)*fmax
-        dP_dtheta3 = np.zeros(I_stokes.shape)*fmax
-        dP_dtheta1[mask] = 1./(I_stokes[mask]*np.sqrt(Q_stokes[mask]**2 + U_stokes[mask]**2))*(Q_stokes[mask]*dQ_dtheta1[mask] + U_stokes[mask]*dQ_dtheta1[mask]) - 1./I_stokes[mask]**2 * dI_dtheta1[mask]
-        dP_dtheta2[mask] = 1./(I_stokes[mask]*np.sqrt(Q_stokes[mask]**2 + U_stokes[mask]**2))*(Q_stokes[mask]*dQ_dtheta2[mask] + U_stokes[mask]*dQ_dtheta2[mask]) - 1./I_stokes[mask]**2 * dI_dtheta2[mask]
-        dP_dtheta3[mask] = 1./(I_stokes[mask]*np.sqrt(Q_stokes[mask]**2 + U_stokes[mask]**2))*(Q_stokes[mask]*dQ_dtheta3[mask] + U_stokes[mask]*dQ_dtheta3[mask]) - 1./I_stokes[mask]**2 * dI_dtheta3[mask]
-        dP_dtheta = np.array([dP_dtheta1, dP_dtheta2, dP_dtheta3])
-
-        dPA_dtheta1 = np.zeros(I_stokes.shape)*fmax
-        dPA_dtheta2 = np.zeros(I_stokes.shape)*fmax
-        dPA_dtheta3 = np.zeros(I_stokes.shape)*fmax
-        dPA_dtheta1[mask] = 90./np.pi*1./(U_stokes[mask]**2 + Q_stokes[mask]**2)*(Q_stokes[mask]*dU_dtheta1[mask] - U_stokes[mask]*dQ_dtheta1[mask])
-        dPA_dtheta2[mask] = 90./np.pi*1./(U_stokes[mask]**2 + Q_stokes[mask]**2)*(Q_stokes[mask]*dU_dtheta2[mask] - U_stokes[mask]*dQ_dtheta2[mask])
-        dPA_dtheta3[mask] = 90./np.pi*1./(U_stokes[mask]**2 + Q_stokes[mask]**2)*(Q_stokes[mask]*dU_dtheta3[mask] - U_stokes[mask]*dQ_dtheta3[mask])
-        dPA_dtheta = np.array([dPA_dtheta1, dPA_dtheta2, dPA_dtheta3])
-
         # Compute the uncertainty associated with the polarizers' orientation (see Kishimoto 1999)
         s_I2_axis = np.sum([dI_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))],axis=0)
         s_Q2_axis = np.sum([dQ_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))],axis=0)
@@ -1196,10 +1167,10 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
             I_stokes, Q_stokes, U_stokes = Stokes_array
             Stokes_cov[0,0], Stokes_cov[1,1], Stokes_cov[2,2] = Stokes_error**2
 
-    return I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta
+    return I_stokes, Q_stokes, U_stokes, Stokes_cov
 
 
-def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta, headers):
+def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers):
     """
     Compute the polarization degree (in %) and angle (in deg) and their
     respective errors from given Stokes parameters.
@@ -1260,14 +1231,6 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta,
     # Propagate previously computed errors
     s_P[mask] = (1/I_stokes[mask])*np.sqrt((Q_stokes[mask]**2*Stokes_cov[1,1][mask] + U_stokes[mask]**2*Stokes_cov[2,2][mask] + 2.*Q_stokes[mask]*U_stokes[mask]*Stokes_cov[1,2][mask])/(Q_stokes[mask]**2 + U_stokes[mask]**2) + ((Q_stokes[mask]/I_stokes[mask])**2 + (U_stokes[mask]/I_stokes[mask])**2)*Stokes_cov[0,0][mask] - 2.*(Q_stokes[mask]/I_stokes[mask])*Stokes_cov[0,1][mask] - 2.*(U_stokes[mask]/I_stokes[mask])*Stokes_cov[0,2][mask])
     s_PA[mask] = (90./(np.pi*(Q_stokes[mask]**2 + U_stokes[mask]**2)))*np.sqrt(U_stokes[mask]**2*Stokes_cov[1,1][mask] + Q_stokes[mask]**2*Stokes_cov[2,2][mask] - 2.*Q_stokes[mask]*U_stokes[mask]*Stokes_cov[1,2][mask])
-    # Compute the uncertainty associated with the polarizers' orientation (see Kishimoto 1999)
-    s_P_axis = np.sqrt(np.sum([dP_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))], axis=0))
-    s_PA_axis = np.sqrt(np.sum([dPA_dtheta[i]**2 * sigma_theta[i]**2 for i in range(len(sigma_theta))], axis=0))
-    # Add axis uncertainties quadratically
-    with warnings.catch_warnings(record=True) as w:
-        s_P[mask] = np.sqrt(s_P[mask]**2 + s_P_axis[mask]**2)
-        s_PA[mask] = np.sqrt(s_PA[mask]**2 + s_PA_axis[mask]**2)
-
     s_P[np.isnan(s_P)] = fmax
     s_PA[np.isnan(s_PA)] = fmax
 
@@ -1303,7 +1266,7 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta,
     return P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P
 
 
-def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dtheta, data_mask, headers, ang, SNRi_cut=None):
+def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, ang, SNRi_cut=None):
     """
     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
@@ -1386,13 +1349,6 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dthet
         new_Stokes_cov[i,i] = np.abs(new_Stokes_cov[i,i])
     center = np.array(new_I_stokes.shape)/2
 
-    #Compute new derivatives of P, PA on rotated parameters
-    new_dP_dtheta = deepcopy(dP_dtheta)
-    new_dPA_dtheta = deepcopy(dPA_dtheta)
-    for i in range(len(dP_dtheta)):
-        new_dP_dtheta[i] = sc_rotate(dP_dtheta[i], ang, reshape=False, cval=0.)
-        new_dPA_dtheta[i] = sc_rotate(dPA_dtheta[i], ang, reshape=False, cval=0.)
-    
     #Update headers to new angle
     new_headers = []
     mrot = np.array([[np.cos(-alpha), -np.sin(-alpha)],
@@ -1435,7 +1391,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, dP_dtheta, dPA_dthet
     new_U_stokes[np.isnan(new_U_stokes)] = 0.
     new_Stokes_cov[np.isnan(new_Stokes_cov)] = fmax
 
-    return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_dP_dtheta, new_dPA_dtheta, new_headers, new_data_mask
+    return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_headers, new_data_mask
 
 def rotate_data(data_array, error_array, data_mask, headers, ang):
     """