correct axis error estimation

src/FOC_reduction.py

@@ -110,8 +110,8 @@ def main():
     # Polarization map output
     figname = 'NGC1068_FOC'         #target/intrument name
     figtype = '_combine_FWHM020_rot'    #additionnal informations
-    SNRp_cut = 20.    #P measurments with SNR>3
-    SNRi_cut = 200   #I measurments with SNR>30, which implies an uncertainty in P of 4.7%.
+    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
 
     ##### Pipeline start

src/lib/Derivative.cpp

@@ -0,0 +1,171 @@
+// Compute the error from the uncertainties in the direction of polarizer's axes
+
+#include <iostream>
+#include <math.h>
+
+
+using namespace std;
+
+const double PI=3.14159265359;	
+
+const double I = 8.42755038425555;
+const double Q = -0.7409864902131016;
+const double U = -1.079689321440609;
+
+const double k1 = 0.986;			//from Kishimoto (1999)
+const double k2 = 0.976;			//from Kishimoto (1999)
+const double k3 = 0.973;			//from Kishimoto (1999)
+
+const double Theta1 = 180*PI/180.;		//radians
+const double Theta2 = 60*PI/180.;		//radians
+const double Theta3 = 120*PI/180.;		//radians
+
+const double sigma_theta_1=3*PI/180.;		//radians
+const double sigma_theta_2=3*PI/180.;		//radians
+const double sigma_theta_3=3*PI/180.;		//radians
+
+
+
+double A()		
+{
+ 	return k2*k3*sin(-2*Theta2+2*Theta3) + k3*k1*sin(-2*Theta3+2*Theta1) + k1*k2*sin(-2*Theta1+2*Theta2);
+}
+
+
+double Ap(int i)		//A'_i
+{
+	if(i==1) return 2*k3*k1*cos(-2*Theta3+2*Theta1) - 2*k1*k2*cos(-2*Theta1+2*Theta2);
+	if(i==2) return -2*k2*k3*cos(-2*Theta2+2*Theta3) + 2*k1*k2*cos(-2*Theta1+2*Theta2);
+	if(i==3) return 2*k2*k3*cos(-2*Theta2+2*Theta3) - 2*k3*k1*cos(-2*Theta3+2*Theta1);
+	else return 0;
+}
+
+
+double Flux(int i)		//f'_i
+{
+	if(i==1) return I-(k1*cos(2*Theta1))*Q-(k1*sin(2*Theta1))*U;
+	if(i==2) return I-(k2*cos(2*Theta2))*Q-(k2*sin(2*Theta2))*U;
+	if(i==3) return I-(k3*cos(2*Theta3))*Q-(k3*sin(2*Theta3))*U;
+	else return 0;
+}
+
+
+double coef(int i, int j)	//a_ij
+{
+	if(i==1 && j==1) return (1./A())*(k2*k3*sin(-2*Theta2+2*Theta3));
+	if(i==1 && j==2) return (1./A())*(k3*k1*sin(-2*Theta3+2*Theta1));
+	if(i==1 && j==3) return (1./A())*(k1*k2*sin(-2*Theta1+2*Theta2));
+	if(i==2 && j==1) return (1./A())*(-k2*sin(2*Theta2)+k3*sin(2*Theta3));
+	if(i==2 && j==2) return (1./A())*(-k3*sin(2*Theta3)+k1*sin(2*Theta1));
+	if(i==2 && j==3) return (1./A())*(-k1*sin(2*Theta1)+k2*sin(2*Theta2));
+	if(i==3 && j==1) return (1./A())*(k2*cos(2*Theta2)-k3*cos(2*Theta3));
+	if(i==3 && j==2) return (1./A())*(k3*cos(2*Theta3)-k1*cos(2*Theta1));
+	if(i==3 && j==3) return (1./A())*(k1*cos(2*Theta1)-k2*cos(2*Theta2));
+	else return 0;	
+}
+
+
+double g()		// I_2 == Q
+{
+	double a,b,c;
+
+	a=coef(2,1)*Flux(1);
+	b=coef(2,2)*Flux(2);
+	c=coef(2,3)*Flux(3);
+	
+ 	return a+b+c;
+}
+
+
+double g_p()	// I'_2 == Q'
+{
+	double f=coef(2,1)*Flux(1)+coef(2,2)*Flux(2)+coef(2,3)*Flux(3);
+	double fp=2*k2*cos(2*Theta2)*(k1*Q*cos(2*Theta1)+k1*U*sin(2*Theta1)-I)+(k1*sin(2*Theta1)-k3*sin(2*Theta3))*(2*k2*Q*sin(2*Theta2)-2*k2*U*cos(2*Theta2))-2*k2*cos(2*Theta2)*(k3*Q*cos(2*Theta3)+k3*U*sin(2*Theta3)-I);
+	return (A()+fp-f*Ap(2))/(pow(A(),2));
+}
+
+double h()		// I_3 == U
+{
+	double a,b,c;
+
+	a=coef(3,1)*Flux(1);
+	b=coef(3,2)*Flux(2);
+	c=coef(3,3)*Flux(3);
+	
+ 	return a+b+c;
+}
+
+
+double h_p()	// I'_3 == U'
+{
+	double f=coef(3,1)*Flux(1)+coef(3,2)*Flux(2)+coef(3,3)*Flux(3);
+	double fp=-2*k3*sin(2*Theta3)*(k1*Q*cos(2*Theta1)+k1*U*sin(2*Theta1)-I)+2*k3*sin(2*Theta3)*(k2*Q*cos(2*Theta2)+k2*U*sin(2*Theta2)-I)+2*k3*(k2*cos(2*Theta2)-k1*cos(2*Theta1))*(U*cos(2*Theta3)-Q*sin(2*Theta3));
+	return (A()+fp-f*Ap(3))/(pow(A(),2));
+}
+
+
+double k()		// I_1 == I
+{
+	double a,b,c;
+
+	a=coef(1,1)*Flux(1);
+	b=coef(1,2)*Flux(2);
+	c=coef(1,3)*Flux(3);
+	
+ 	return a+b+c;
+}
+
+
+double k_p()	// I'_1 == I'
+{
+	double f=coef(1,1)*Flux(1)+coef(1,2)*Flux(2)+coef(1,3)*Flux(3);
+	double fp=2*k1*k2*k3*sin(2*Theta2-2*Theta3)*(U*cos(2*Theta1)-Q*sin(2*Theta1))-2*k1*k3*cos(2*Theta1-2*Theta3)*(k2*Q*cos(2*Theta2)+k2*U*sin(2*Theta2)-I)+2*k1*k2*cos(2*Theta1-2*Theta2)*(k3*Q*cos(2*Theta3)+k3*U*sin(2*Theta3)-I);
+	return (A()+fp-f*Ap(1))/(pow(A(),2));
+}
+
+
+double dTheta_i(int i)		
+{
+	double a,b,c;
+
+	if(i==1)			//partial derivative with i=1
+	{
+		a = g()*g_p();
+		b = k();
+		c = sqrt(pow(g(),2)+pow(h(),2));
+		return a/(b*c);
+	}
+	if(i==2)			//partial derivative with i=2
+	{
+		a = h()*h_p();
+		b = k();
+		c = sqrt(pow(g(),2)+pow(h(),2));
+		return a/(b*c);
+	}
+	if(i==3)			//partial derivative with i=3
+	{
+		a = k_p();
+		b = sqrt(pow(g(),2)+pow(h(),2));
+		c = pow(k(),2);
+		return -a*b/c;
+	}
+	else return 0;
+}
+
+
+
+int main()
+{
+	double sigma_stat_P;
+	sigma_stat_P = 0;
+
+	sigma_stat_P=pow(dTheta_i(1),2)*pow(sigma_theta_1,2)+pow(dTheta_i(2),2)*pow(sigma_theta_2,2)+pow(dTheta_i(3),2)*pow(sigma_theta_3,2);
+
+	cout << "P: " << sqrt(Q*Q+U*U)/I << " +/- " << sqrt(sigma_stat_P) << endl;
+
+	return 0;
+}
+
+
+
+

src/lib/plots.py

@@ -359,7 +359,7 @@ def polarization_map(Stokes, rectangle=None, SNRp_cut=3., SNRi_cut=30., step_vec
     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_err = P_diluted_err/(2.*P_diluted)*180./np.pi
 
-    ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(pivot_wav,sci_not(I_diluted*convert_flux,I_diluted_err*convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_diluted*100.,P_diluted_err*100.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_diluted,PA_diluted_err), color='white', fontsize=16, xy=(0.01, 0.92), xycoords='axes fraction')
+    ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(pivot_wav,sci_not(I_int*convert_flux,I_int_err*convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_int*100.,P_int_err*100.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_int,PA_int_err), color='white', fontsize=16, xy=(0.01, 0.92), xycoords='axes fraction')
 
     ax.coords.grid(True, color='white', ls='dotted', alpha=0.5)
     ax.coords[0].set_axislabel('Right Ascension (J2000)')

src/lib/reduction.py

@@ -1106,21 +1106,22 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
         pol_eff[1] = pol_efficiency['pol60']
         pol_eff[2] = pol_efficiency['pol120']
 
-        theta = np.array([np.pi, np.pi/3., 2.*np.pi/3.])
+        theta = np.array([180.*np.pi/180., 60.*np.pi/180., 120.*np.pi/180.])
+        sigma_theta = np.array([3.*np.pi/180., 3.*np.pi/180., 3.*np.pi/180.])
         pol_flux = 2.*np.array([pol0/transmit[0], pol60/transmit[1], pol120/transmit[2]])
 
-        norm = pol_eff[1]*pol_eff[2]*np.sin(-2.*theta[1]+2.*theta[2]) \
+        A = pol_eff[1]*pol_eff[2]*np.sin(-2.*theta[1]+2.*theta[2]) \
                 + pol_eff[2]*pol_eff[0]*np.sin(-2.*theta[2]+2.*theta[0]) \
                 + pol_eff[0]*pol_eff[1]*np.sin(-2.*theta[0]+2.*theta[1])
-        a_stokes = np.zeros((3,3))
+        coeff_stokes = np.zeros((3,3))
         for i in range(3):
-            a_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])/norm
-            a_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]))/norm
-            a_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]))/norm
+            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])/A
+            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]))/A
+            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]))/A
 
-        I_stokes = np.sum([a_stokes[0,i]*pol_flux[i] for i in range(3)], axis=0)
-        Q_stokes = np.sum([a_stokes[1,i]*pol_flux[i] for i in range(3)], axis=0)
-        U_stokes = np.sum([a_stokes[2,i]*pol_flux[i] for i in range(3)], axis=0)
+        I_stokes = np.sum([coeff_stokes[0,i]*pol_flux[i] for i in range(3)], axis=0)
+        Q_stokes = np.sum([coeff_stokes[1,i]*pol_flux[i] for i in range(3)], axis=0)
+        U_stokes = np.sum([coeff_stokes[2,i]*pol_flux[i] for i in range(3)], axis=0)
 
         # Remove nan
         I_stokes[np.isnan(I_stokes)]=0.
@@ -1135,27 +1136,35 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
 
         #Stokes covariance matrix
         Stokes_cov = np.zeros((3,3,I_stokes.shape[0],I_stokes.shape[1]))
-        Stokes_cov[0,0] = a_stokes[0,0]**2*pol_cov[0,0]+a_stokes[0,1]**2*pol_cov[1,1]+a_stokes[0,2]**2*pol_cov[2,2] + 2*(a_stokes[0,0]*a_stokes[0,1]*pol_cov[0,1]+a_stokes[0,0]*a_stokes[0,2]*pol_cov[0,2]+a_stokes[0,1]*a_stokes[0,2]*pol_cov[1,2])
-        Stokes_cov[1,1] = a_stokes[1,0]**2*pol_cov[0,0]+a_stokes[1,1]**2*pol_cov[1,1]+a_stokes[1,2]**2*pol_cov[2,2] + 2*(a_stokes[1,0]*a_stokes[1,1]*pol_cov[0,1]+a_stokes[1,0]*a_stokes[1,2]*pol_cov[0,2]+a_stokes[1,1]*a_stokes[1,2]*pol_cov[1,2])
-        Stokes_cov[2,2] = a_stokes[2,0]**2*pol_cov[0,0]+a_stokes[2,1]**2*pol_cov[1,1]+a_stokes[2,2]**2*pol_cov[2,2] + 2*(a_stokes[2,0]*a_stokes[2,1]*pol_cov[0,1]+a_stokes[2,0]*a_stokes[2,2]*pol_cov[0,2]+a_stokes[2,1]*a_stokes[2,2]*pol_cov[1,2])
-        Stokes_cov[0,1] = Stokes_cov[1,0] = a_stokes[0,0]*a_stokes[1,0]*pol_cov[0,0]+a_stokes[0,1]*a_stokes[1,1]*pol_cov[1,1]+a_stokes[0,2]*a_stokes[1,2]*pol_cov[2,2]+(a_stokes[0,0]*a_stokes[1,1]+a_stokes[1,0]*a_stokes[0,1])*pol_cov[0,1]+(a_stokes[0,0]*a_stokes[1,2]+a_stokes[1,0]*a_stokes[0,2])*pol_cov[0,2]+(a_stokes[0,1]*a_stokes[1,2]+a_stokes[1,1]*a_stokes[0,2])*pol_cov[1,2]
-        Stokes_cov[0,2] = Stokes_cov[2,0] = a_stokes[0,0]*a_stokes[2,0]*pol_cov[0,0]+a_stokes[0,1]*a_stokes[2,1]*pol_cov[1,1]+a_stokes[0,2]*a_stokes[2,2]*pol_cov[2,2]+(a_stokes[0,0]*a_stokes[2,1]+a_stokes[2,0]*a_stokes[0,1])*pol_cov[0,1]+(a_stokes[0,0]*a_stokes[2,2]+a_stokes[2,0]*a_stokes[0,2])*pol_cov[0,2]+(a_stokes[0,1]*a_stokes[2,2]+a_stokes[2,1]*a_stokes[0,2])*pol_cov[1,2]
-        Stokes_cov[1,2] = Stokes_cov[2,1] = a_stokes[1,0]*a_stokes[2,0]*pol_cov[0,0]+a_stokes[1,1]*a_stokes[2,1]*pol_cov[1,1]+a_stokes[1,2]*a_stokes[2,2]*pol_cov[2,2]+(a_stokes[1,0]*a_stokes[2,1]+a_stokes[2,0]*a_stokes[1,1])*pol_cov[0,1]+(a_stokes[1,0]*a_stokes[2,2]+a_stokes[2,0]*a_stokes[1,2])*pol_cov[0,2]+(a_stokes[1,1]*a_stokes[2,2]+a_stokes[2,1]*a_stokes[1,2])*pol_cov[1,2]
+        Stokes_cov[0,0] = coeff_stokes[0,0]**2*pol_cov[0,0]+coeff_stokes[0,1]**2*pol_cov[1,1]+coeff_stokes[0,2]**2*pol_cov[2,2] + 2*(coeff_stokes[0,0]*coeff_stokes[0,1]*pol_cov[0,1]+coeff_stokes[0,0]*coeff_stokes[0,2]*pol_cov[0,2]+coeff_stokes[0,1]*coeff_stokes[0,2]*pol_cov[1,2])
+        Stokes_cov[1,1] = coeff_stokes[1,0]**2*pol_cov[0,0]+coeff_stokes[1,1]**2*pol_cov[1,1]+coeff_stokes[1,2]**2*pol_cov[2,2] + 2*(coeff_stokes[1,0]*coeff_stokes[1,1]*pol_cov[0,1]+coeff_stokes[1,0]*coeff_stokes[1,2]*pol_cov[0,2]+coeff_stokes[1,1]*coeff_stokes[1,2]*pol_cov[1,2])
+        Stokes_cov[2,2] = coeff_stokes[2,0]**2*pol_cov[0,0]+coeff_stokes[2,1]**2*pol_cov[1,1]+coeff_stokes[2,2]**2*pol_cov[2,2] + 2*(coeff_stokes[2,0]*coeff_stokes[2,1]*pol_cov[0,1]+coeff_stokes[2,0]*coeff_stokes[2,2]*pol_cov[0,2]+coeff_stokes[2,1]*coeff_stokes[2,2]*pol_cov[1,2])
+        Stokes_cov[0,1] = Stokes_cov[1,0] = coeff_stokes[0,0]*coeff_stokes[1,0]*pol_cov[0,0]+coeff_stokes[0,1]*coeff_stokes[1,1]*pol_cov[1,1]+coeff_stokes[0,2]*coeff_stokes[1,2]*pol_cov[2,2]+(coeff_stokes[0,0]*coeff_stokes[1,1]+coeff_stokes[1,0]*coeff_stokes[0,1])*pol_cov[0,1]+(coeff_stokes[0,0]*coeff_stokes[1,2]+coeff_stokes[1,0]*coeff_stokes[0,2])*pol_cov[0,2]+(coeff_stokes[0,1]*coeff_stokes[1,2]+coeff_stokes[1,1]*coeff_stokes[0,2])*pol_cov[1,2]
+        Stokes_cov[0,2] = Stokes_cov[2,0] = coeff_stokes[0,0]*coeff_stokes[2,0]*pol_cov[0,0]+coeff_stokes[0,1]*coeff_stokes[2,1]*pol_cov[1,1]+coeff_stokes[0,2]*coeff_stokes[2,2]*pol_cov[2,2]+(coeff_stokes[0,0]*coeff_stokes[2,1]+coeff_stokes[2,0]*coeff_stokes[0,1])*pol_cov[0,1]+(coeff_stokes[0,0]*coeff_stokes[2,2]+coeff_stokes[2,0]*coeff_stokes[0,2])*pol_cov[0,2]+(coeff_stokes[0,1]*coeff_stokes[2,2]+coeff_stokes[2,1]*coeff_stokes[0,2])*pol_cov[1,2]
+        Stokes_cov[1,2] = Stokes_cov[2,1] = coeff_stokes[1,0]*coeff_stokes[2,0]*pol_cov[0,0]+coeff_stokes[1,1]*coeff_stokes[2,1]*pol_cov[1,1]+coeff_stokes[1,2]*coeff_stokes[2,2]*pol_cov[2,2]+(coeff_stokes[1,0]*coeff_stokes[2,1]+coeff_stokes[2,0]*coeff_stokes[1,1])*pol_cov[0,1]+(coeff_stokes[1,0]*coeff_stokes[2,2]+coeff_stokes[2,0]*coeff_stokes[1,2])*pol_cov[0,2]+(coeff_stokes[1,1]*coeff_stokes[2,2]+coeff_stokes[2,1]*coeff_stokes[1,2])*pol_cov[1,2]
 
-        C1 = 2.*pol_eff[0]*pol_eff[1]*pol_eff[2]/norm
-        dI_dtheta1 = C1*(np.cos(-2.*theta[2]+2.*theta[0])/pol_eff[1]*(pol_flux[1]-I_stokes) - np.cos(-2.*theta[0]+2.*theta[1])/pol_eff[2]*(pol_flux[2]-I_stokes))
-        dI_dtheta2 = C1*(np.cos(-2.*theta[0]+2.*theta[1])/pol_eff[2]*(pol_flux[2]-I_stokes) - np.cos(-2.*theta[1]+2.*theta[2])/pol_eff[0]*(pol_flux[0]-I_stokes))
-        dI_dtheta3 = C1*(np.cos(-2.*theta[1]+2.*theta[2])/pol_eff[0]*(pol_flux[0]-I_stokes) - np.cos(-2.*theta[2]+2.*theta[0])/pol_eff[1]*(pol_flux[1]-I_stokes))
-        dQ_dtheta1 = C1*((np.cos(2.*theta[0])*pol_flux[1]-np.cos(2.*theta[0])*pol_flux[2])/(pol_eff[1]*pol_eff[2]) - (np.cos(-2.*theta[2]+2.*theta[0])/pol_eff[1]-np.cos(-2.*theta[0]+2.*theta[1])/pol_eff[2])*Q_stokes)
-        dQ_dtheta2 = C1*((np.cos(2.*theta[1])*pol_flux[2]-np.cos(2.*theta[1])*pol_flux[0])/(pol_eff[0]*pol_eff[2]) - (np.cos(-2.*theta[0]+2.*theta[1])/pol_eff[2]-np.cos(-2.*theta[1]+2.*theta[2])/pol_eff[0])*Q_stokes)
-        dQ_dtheta3 = C1*((np.cos(2.*theta[2])*pol_flux[0]-np.cos(2.*theta[2])*pol_flux[1])/(pol_eff[0]*pol_eff[1]) - (np.cos(-2.*theta[1]+2.*theta[2])/pol_eff[0]-np.cos(-2.*theta[2]+2.*theta[0])/pol_eff[1])*Q_stokes)
-        dU_dtheta1 = C1*((np.sin(2.*theta[0])*pol_flux[1]-np.sin(2.*theta[1])*pol_flux[2])/(pol_eff[1]*pol_eff[2]) - (np.cos(-2.*theta[2]+2.*theta[0])/pol_eff[1]-np.cos(-2.*theta[0]+2.*theta[1])/pol_eff[2])*U_stokes)
-        dU_dtheta2 = C1*((np.sin(2.*theta[1])*pol_flux[2]-np.sin(2.*theta[1])*pol_flux[0])/(pol_eff[0]*pol_eff[2]) - (np.cos(-2.*theta[0]+2.*theta[1])/pol_eff[2]-np.cos(-2.*theta[1]+2.*theta[2])/pol_eff[0])*U_stokes)
-        dU_dtheta3 = C1*((np.sin(2.*theta[2])*pol_flux[0]-np.sin(2.*theta[2])*pol_flux[1])/(pol_eff[0]*pol_eff[1]) - (np.cos(-2.*theta[1]+2.*theta[2])/pol_eff[0]-np.cos(-2.*theta[2]+2.*theta[0])/pol_eff[1])*U_stokes)
+        dI_dtheta1 = 2.*pol_eff[0]/A*(pol_eff[2]*np.cos(-2.*theta[2]+2.*theta[0])*(pol_flux[1]-I_stokes) - pol_eff[1]*np.cos(-2.*theta[0]+2.*theta[1])*(pol_flux[2]-I_stokes) + A*coeff_stokes[0,0]*(np.sin(2.*theta[0]*Q_stokes) - np.cos(2.*theta[0]*U_stokes)))
+        dI_dtheta2 = 2.*pol_eff[1]/A*(pol_eff[0]*np.cos(-2.*theta[0]+2.*theta[1])*(pol_flux[2]-I_stokes) - pol_eff[2]*np.cos(-2.*theta[1]+2.*theta[2])*(pol_flux[0]-I_stokes) + A*coeff_stokes[0,1]*(np.sin(2.*theta[1]*Q_stokes) - np.cos(2.*theta[1]*U_stokes)))
+        dI_dtheta3 = 2.*pol_eff[2]/A*(pol_eff[1]*np.cos(-2.*theta[1]+2.*theta[2])*(pol_flux[0]-I_stokes) - pol_eff[0]*np.cos(-2.*theta[2]+2.*theta[0])*(pol_flux[1]-I_stokes) + A*coeff_stokes[0,2]*(np.sin(2.*theta[2]*Q_stokes) - np.cos(2.*theta[2]*U_stokes)))
+        dQ_dtheta1 = 2.*pol_eff[0]/A*(np.cos(2.*theta[0])*(pol_flux[1]-pol_flux[2]) - (pol_eff[2]*np.cos(-2.*theta[2]+2.*theta[0]) - pol_eff[1]*np.cos(-2.*theta[0]+2.*theta[1]))*Q_stokes + A*coeff_stokes[1,0]*(np.sin(2.*theta[0]*Q_stokes) - np.cos(2.*theta[0]*U_stokes)))
+        dQ_dtheta2 = 2.*pol_eff[1]/A*(np.cos(2.*theta[1])*(pol_flux[2]-pol_flux[0]) - (pol_eff[0]*np.cos(-2.*theta[0]+2.*theta[1]) - pol_eff[2]*np.cos(-2.*theta[1]+2.*theta[2]))*Q_stokes + A*coeff_stokes[1,1]*(np.sin(2.*theta[1]*Q_stokes) - np.cos(2.*theta[1]*U_stokes)))
+        dQ_dtheta3 = 2.*pol_eff[2]/A*(np.cos(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]))*Q_stokes + A*coeff_stokes[1,2]*(np.sin(2.*theta[2]*Q_stokes) - np.cos(2.*theta[2]*U_stokes)))
+        dU_dtheta1 = 2.*pol_eff[0]/A*(np.sin(2.*theta[0])*(pol_flux[1]-pol_flux[2]) - (pol_eff[2]*np.cos(-2.*theta[2]+2.*theta[0]) - pol_eff[1]*np.cos(-2.*theta[0]+2.*theta[1]))*U_stokes + A*coeff_stokes[2,0]*(np.sin(2.*theta[0]*Q_stokes) - np.cos(2.*theta[0]*U_stokes)))
+        dU_dtheta2 = 2.*pol_eff[1]/A*(np.sin(2.*theta[1])*(pol_flux[2]-pol_flux[0]) - (pol_eff[0]*np.cos(-2.*theta[0]+2.*theta[1]) - pol_eff[2]*np.cos(-2.*theta[1]+2.*theta[2]))*U_stokes + A*coeff_stokes[2,1]*(np.sin(2.*theta[1]*Q_stokes) - np.cos(2.*theta[1]*U_stokes)))
+        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 + A*coeff_stokes[2,2]*(np.sin(2.*theta[2]*Q_stokes) - np.cos(2.*theta[2]*U_stokes)))
 
-        #Stokes_cov[0,0] += (dI_dtheta1 + dI_dtheta2 + dI_dtheta3)**2*3.*np.pi/180.
-        #Stokes_cov[1,1] += (dQ_dtheta1 + dQ_dtheta2 + dQ_dtheta3)**2*3.*np.pi/180.
-        #Stokes_cov[2,2] += (dU_dtheta1 + dU_dtheta2 + dU_dtheta3)**2*3.*np.pi/180.
+        Stokes_cov[0,0] += (dI_dtheta1**2*sigma_theta[0]**2 + dI_dtheta2**2*sigma_theta[1]**2 + dI_dtheta3**2*sigma_theta[2]**2)
+        Stokes_cov[1,1] += (dQ_dtheta1**2*sigma_theta[0]**2 + dQ_dtheta2**2*sigma_theta[1]**2 + dQ_dtheta3**2*sigma_theta[2]**2)
+        Stokes_cov[2,2] += (dU_dtheta1**2*sigma_theta[0]**2 + dU_dtheta2**2*sigma_theta[1]**2 + dU_dtheta3**2*sigma_theta[2]**2)
+        plt.imshow(np.abs(Stokes_cov[0,0]/I_stokes)*100., origin='lower')
+        plt.colorbar()
+        plt.show()
+        plt.imshow(np.abs(Stokes_cov[1,1]/Q_stokes)*100., origin='lower')
+        plt.colorbar()
+        plt.show()
+        plt.imshow(np.abs(Stokes_cov[2,2]/U_stokes)*100., origin='lower')
+        plt.colorbar()
+        plt.show()
 
         if not(FWHM is None) and (smoothing.lower() in ['gaussian_after','gauss_after']):
             Stokes_array = np.array([I_stokes, Q_stokes, U_stokes])