diff --git a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot.png b/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot.png deleted file mode 100644 index 76e609c..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_IQU.png b/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_IQU.png deleted file mode 100644 index eee2321..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_IQU.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_P.png b/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_P.png deleted file mode 100644 index 438498b..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_P.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_P_err.png deleted file mode 100644 index ead8731..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_P_err.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_P_flux.png deleted file mode 100644 index 0a52bf9..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_P_flux.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_SNRi.png deleted file mode 100644 index c780e59..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_SNRi.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_SNRp.png deleted file mode 100644 index f6c8309..0000000 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_2_combine_FWHM020_rot_SNRp.png and /dev/null differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot.png index cd7b911..0c7e32c 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P.png index 5082534..9305f07 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P_err.png index 004e35c..69e0bb0 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P_err.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P_flux.png index cfd5da8b..0be4f8f 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P_flux.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_P_flux.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRi.png index aa5b637..5d4cec2 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRi.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRp.png index 639786c..d5214f1 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRp.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM010_rot_SNRp.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot.png index 84ae9e2..0c34ef5 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P.png index 0c1b977..98b5b4c 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_err.png index 5b19875..732fc19 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_err.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_flux.png index e8b5320..da1d12e 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_flux.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_flux.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRi.png index cee85cf..be899b2 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRi.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRp.png index ebd3494..c72c900 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRp.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRp.png differ diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 57a7239..28f2433 100755 --- a/src/FOC_reduction.py +++ b/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 diff --git a/src/lib/Derivative.cpp b/src/lib/Derivative.cpp new file mode 100644 index 0000000..10a290c --- /dev/null +++ b/src/lib/Derivative.cpp @@ -0,0 +1,171 @@ +// Compute the error from the uncertainties in the direction of polarizer's axes + +#include +#include + + +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; +} + + + + diff --git a/src/lib/plots.py b/src/lib/plots.py index e5d5060..4dc2bad 100755 --- a/src/lib/plots.py +++ b/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)') diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 344390c..751db03 100755 --- a/src/lib/reduction.py +++ b/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])