diff --git a/FOC_Capetti_test_NGC1068/FOC_Capetti_test_NGC1068.rar b/FOC_Capetti_test_NGC1068/FOC_Capetti_test_NGC1068.rar new file mode 100755 index 0000000..f9b9f44 Binary files /dev/null and b/FOC_Capetti_test_NGC1068/FOC_Capetti_test_NGC1068.rar differ diff --git a/FOC_Capetti_test_NGC1068/FOC_reduction.py b/FOC_Capetti_test_NGC1068/FOC_reduction.py new file mode 100755 index 0000000..06dd2b2 --- /dev/null +++ b/FOC_Capetti_test_NGC1068/FOC_reduction.py @@ -0,0 +1,256 @@ +from pylab import * +import numpy as np +import matplotlib.pyplot as plt +from astropy.io import fits +from astropy.wcs import WCS +from aplpy import FITSFigure +import scipy.ndimage +import os as os + +plt.close('all') + +def bin_ndarray(ndarray, new_shape, operation='sum'): + """ + Bins an ndarray in all axes based on the target shape, by summing or + averaging. + + Number of output dimensions must match number of input dimensions. + + Example + ------- + >>> m = np.arange(0,100,1).reshape((10,10)) + >>> n = bin_ndarray(m, new_shape=(5,5), operation='sum') + >>> print(n) + + [[ 22 30 38 46 54] + [102 110 118 126 134] + [182 190 198 206 214] + [262 270 278 286 294] + [342 350 358 366 374]] + + """ + if not operation.lower() in ['sum', 'mean', 'average', 'avg']: + raise ValueError("Operation not supported.") + if ndarray.ndim != len(new_shape): + raise ValueError("Shape mismatch: {} -> {}".format(ndarray.shape, + new_shape)) + compression_pairs = [(d, c//d) for d,c in zip(new_shape, + ndarray.shape)] + flattened = [l for p in compression_pairs for l in p] + ndarray = ndarray.reshape(flattened) + for i in range(len(new_shape)): + if operation.lower() == "sum": + ndarray = ndarray.sum(-1*(i+1)) + elif operation.lower() in ["mean", "average", "avg"]: + ndarray = ndarray.mean(-1*(i+1)) + return ndarray + +def plots(ax,data,vmax,vmin): + ax.imshow(data,vmax=vmax,vmin=vmin,origin='lower') + + + +### User input + +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'] + +#Centroids +#The centroids should be estimated by cross-correlating the images. +#Here I used the position of the central source for each file as the reference pixel position. +ycen_0 = [304,306,303,296,295,295,294,305,304] +xcen_0 = [273,274,273,276,274,274,274,272,271] +#size, in pixels, of the FOV centered in x_cen_0,y_cen_0 +Dx = 200 +Dy = 200 + +#set the new image size (Dxy x Dxy pixels binning) +Dxy = 5 +new_shape = (Dx//Dxy,Dy//Dxy) + +#figures +#test alignment +vmin = 0 +vmax = 6 +font_size=40.0 +label_size=20. +lw = 3. + +#pol. map +SNRp_cut = 3 #P measumentes with SNR>3 +SNRi_cut = 30 #I measuremntes with SNR>30, which implies an uncerrtainty in P of 4.7%. +scalevec = 0.05 #length of vectors in pol. map +step_vec = 1 #plot all vectors in the array. if step_vec = 2, then every other vector will be plotted +vec_legend = 10 #% pol for legend +figname = 'NGC1068_FOC.png' + + +### SCRIPT ### +### Step 1. Check input images before data reduction +#this step is very simplistic. +#Here I used the position of the central source for each file as the +#reference pixel position. +#The centroids should be estimated by cross-correlating the images, +#and compare with the simplistic approach of using the peak pixel of the +#object as the reference pixel. + + +fig,axs = plt.subplots(3,3,figsize=(30,30),dpi=200,sharex=True,sharey=True) + +for jj, a in enumerate(axs.flatten()): + img = fits.open(infiles[jj]) + ima = img[0].data + ima = ima[ycen_0[jj]-Dy:ycen_0[jj]+Dy,xcen_0[jj]-Dx:xcen_0[jj]+Dx] + ima = bin_ndarray(ima,new_shape=new_shape,operation='sum') #binning + exptime = img[0].header['EXPTIME'] + fil = img[0].header['FILTNAM1'] + ima = ima/exptime + globals()['ima_%s' % jj] = ima + #plots + plots(a,ima,vmax=vmax,vmin=vmin) + #position of centroid + a.plot([ima.shape[1]/2,ima.shape[1]/2],[0,ima.shape[0]-1],lw=1,color='black') + a.plot([0,ima.shape[1]-1],[ima.shape[1]/2,ima.shape[1]/2],lw=1,color='black') + a.text(2,2,infiles[jj][0:8],color='white',fontsize=10) + a.text(2,5,fil,color='white',fontsize=30) + a.text(ima.shape[1]-20,1,exptime,color='white',fontsize=20) +fig.subplots_adjust(hspace=0,wspace=0) +fig.savefig('test_alignment.png',dpi=300) +os.system('open test_alignment.png') + + + +### Step 2. average of all images for a single polarizer to have them in the same units e/s. +pol0 = (ima_0 + ima_1 + ima_2)/3. +pol60 = (ima_3 + ima_4 + ima_5 + ima_6)/4. +pol120 = (ima_7 + ima_8)/2. + +fig1,(ax1,ax2,ax3) = plt.subplots(1,3,figsize=(26,8),dpi=200) +CF = ax1.imshow(pol0,vmin=vmin,vmax=vmax,origin='lower') +cbar = plt.colorbar(CF,ax=ax1) +cbar.ax.tick_params(labelsize=20) +ax1.tick_params(axis='both', which='major', labelsize=20) +ax1.text(2,2,'POL0',color='white',fontsize=10) + +CF = ax2.imshow(pol60,vmin=vmin,vmax=vmax,origin='lower') +cbar = plt.colorbar(CF,ax=ax2) +cbar.ax.tick_params(labelsize=20) +ax2.tick_params(axis='both', which='major', labelsize=20) +ax2.text(2,2,'POL60',color='white',fontsize=10) + +CF = ax3.imshow(pol120,vmin=vmin,vmax=vmax,origin='lower') +cbar = plt.colorbar(CF,ax=ax3) +cbar.ax.tick_params(labelsize=20) +ax3.tick_params(axis='both', which='major', labelsize=20) +ax3.text(2,2,'POL120',color='white',fontsize=10) +fig1.savefig('test_combinePol.png',dpi=300) +os.system('open test_combinePol.png') + + +### Step 3. Compute Stokes IQU, P, PA, PI +#Stokes parameters +I_stokes = (2./3.)*(pol0 + pol60 + pol120) +Q_stokes = (2./3.)*(2*pol0 - pol60 - pol120) +U_stokes = (2./np.sqrt(3.))*(pol60 - pol120) + +#Remove nan +I_stokes[np.isnan(I_stokes)]=0. +Q_stokes[np.isnan(Q_stokes)]=0. +U_stokes[np.isnan(U_stokes)]=0. + +#Polarimetry +PI = np.sqrt(Q_stokes*Q_stokes + U_stokes*U_stokes) +P = PI/I_stokes*100 +PA = 0.5*arctan2(U_stokes,Q_stokes)*180./np.pi+90 +s_P = np.sqrt(2.)*(I_stokes)**(-0.5) +s_PA = s_P/(P/100.)*180./np.pi + +fig2,((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2,3,figsize=(40,20),dpi=200) +CF = ax1.imshow(I_stokes,origin='lower') +cbar = plt.colorbar(CF,ax=ax1) +cbar.ax.tick_params(labelsize=20) +ax1.tick_params(axis='both', which='major', labelsize=20) +ax1.text(2,2,'I',color='white',fontsize=10) + +CF = ax2.imshow(Q_stokes,origin='lower') +cbar = plt.colorbar(CF,ax=ax2) +cbar.ax.tick_params(labelsize=20) +ax2.tick_params(axis='both', which='major', labelsize=20) +ax2.text(2,2,'Q',color='white',fontsize=10) + +CF = ax3.imshow(U_stokes,origin='lower') +cbar = plt.colorbar(CF,ax=ax3) +cbar.ax.tick_params(labelsize=20) +ax3.tick_params(axis='both', which='major', labelsize=20) +ax3.text(2,2,'U',color='white',fontsize=10) + +v = np.linspace(0,40,50) +CF = ax4.imshow(P,origin='lower',vmin=0,vmax=40) +cbar = plt.colorbar(CF,ax=ax4) +cbar.ax.tick_params(labelsize=20) +ax4.tick_params(axis='both', which='major', labelsize=20) +ax4.text(2,2,'P',color='white',fontsize=10) + +CF = ax5.imshow(PA,origin='lower',vmin=0,vmax=180) +cbar = plt.colorbar(CF,ax=ax5) +cbar.ax.tick_params(labelsize=20) +ax5.tick_params(axis='both', which='major', labelsize=20) +ax5.text(2,2,'PA',color='white',fontsize=10) + +CF = ax6.imshow(PI,origin='lower') +cbar = plt.colorbar(CF,ax=ax6) +cbar.ax.tick_params(labelsize=20) +ax6.tick_params(axis='both', which='major', labelsize=20) +ax6.text(2,2,'PI',color='white',fontsize=10) + +fig2.savefig('test_Stokes.png',dpi=300) +os.system('open test_Stokes.png') + +### Step 4. Binning and smoothing +#Images can be binned and smoothed to improve SNR. This step can also be done +#using the PolX images. + + +### Step 5. Roate images to have North up +#Images needs to be reprojected to have North up. +#this procedure implies to rotate the Stokes QU using a rotation matrix + + +### STEP 6. image to FITS with updated WCS +new_wcs = WCS(naxis=2) +new_wcs.wcs.crpix = [I_stokes.shape[0]/2, I_stokes.shape[1]/2] +new_wcs.wcs.crval = [img[0].header['CRVAL1'], img[0].header['CRVAL2']] +new_wcs.wcs.cunit = ["deg", "deg"] +new_wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"] +new_wcs.wcs.cdelt = [img[0].header['CD1_1']*Dxy, img[0].header['CD1_2']*Dxy] + +#hdu_ori = WCS(img[0]) +stkI=fits.PrimaryHDU(data=I_stokes,header=new_wcs.to_header()) +pol=fits.PrimaryHDU(data=P,header=new_wcs.to_header()) +pang=fits.PrimaryHDU(data=PA,header=new_wcs.to_header()) +pol_err=fits.PrimaryHDU(data=s_P,header=new_wcs.to_header()) +pang_err=fits.PrimaryHDU(data=s_PA,header=new_wcs.to_header()) + + +### STEP 7. polarization map +#quality cuts +pxscale = stkI.header['CDELT1'] + +#apply quality cuts +SNRp = pol.data/pol_err.data +pol.data[SNRp < SNRp_cut] = np.nan + +SNRi = stkI.data/np.std(stkI.data[0:10,0:10]) +pol.data[SNRi < SNRi_cut] = np.nan + +fig = plt.figure(figsize=(11,10)) +gc = FITSFigure(stkI,figure=fig) +gc.show_contour(np.log10(SNRi),levels=np.linspace(np.log10(SNRi_cut),np.max(np.log10(SNRi)),20),\ + filled=True,cmap='magma') +gc.show_vectors(pol,pang,scale=scalevec,step=step_vec,color='white',linewidth=1.0) + +fig.savefig(figname,dpi=300) +os.system('open '+figname) + + diff --git a/FOC_Capetti_test_NGC1068/NGC1068_FOC.png b/FOC_Capetti_test_NGC1068/NGC1068_FOC.png new file mode 100755 index 0000000..76da21e Binary files /dev/null and b/FOC_Capetti_test_NGC1068/NGC1068_FOC.png differ diff --git a/FOC_Capetti_test_NGC1068/test_Stokes.png b/FOC_Capetti_test_NGC1068/test_Stokes.png new file mode 100755 index 0000000..ad51433 Binary files /dev/null and b/FOC_Capetti_test_NGC1068/test_Stokes.png differ diff --git a/FOC_Capetti_test_NGC1068/test_alignment.png b/FOC_Capetti_test_NGC1068/test_alignment.png new file mode 100755 index 0000000..ef2d020 Binary files /dev/null and b/FOC_Capetti_test_NGC1068/test_alignment.png differ diff --git a/FOC_Capetti_test_NGC1068/test_combinePol.png b/FOC_Capetti_test_NGC1068/test_combinePol.png new file mode 100755 index 0000000..2938603 Binary files /dev/null and b/FOC_Capetti_test_NGC1068/test_combinePol.png differ diff --git a/cfitsio-3.49.tar.gz b/cfitsio-3.49.tar.gz new file mode 100755 index 0000000..be8bfe3 Binary files /dev/null and b/cfitsio-3.49.tar.gz differ diff --git a/plots/3C405_x136060/3C405_FOC.png b/plots/3C405_x136060/3C405_FOC.png index dbd9734..8207017 100755 Binary files a/plots/3C405_x136060/3C405_FOC.png and b/plots/3C405_x136060/3C405_FOC.png differ diff --git a/plots/3C405_x136060/3C405_FOC_IQU.png b/plots/3C405_x136060/3C405_FOC_IQU.png index 73dc614..645e1ca 100755 Binary files a/plots/3C405_x136060/3C405_FOC_IQU.png and b/plots/3C405_x136060/3C405_FOC_IQU.png differ diff --git a/plots/3C405_x136060/3C405_FOC_P.png b/plots/3C405_x136060/3C405_FOC_P.png index c0b2bb7..413d2a0 100755 Binary files a/plots/3C405_x136060/3C405_FOC_P.png and b/plots/3C405_x136060/3C405_FOC_P.png differ diff --git a/plots/3C405_x136060/3C405_FOC_P_err.png b/plots/3C405_x136060/3C405_FOC_P_err.png index 3103e98..8c5aab5 100755 Binary files a/plots/3C405_x136060/3C405_FOC_P_err.png and b/plots/3C405_x136060/3C405_FOC_P_err.png differ diff --git a/plots/3C405_x136060/3C405_FOC_P_flux.png b/plots/3C405_x136060/3C405_FOC_P_flux.png new file mode 100644 index 0000000..bd7d4d0 Binary files /dev/null and b/plots/3C405_x136060/3C405_FOC_P_flux.png differ diff --git a/plots/3C405_x136060/3C405_FOC_SNRi.png b/plots/3C405_x136060/3C405_FOC_SNRi.png index a159ed1..1aa13f1 100755 Binary files a/plots/3C405_x136060/3C405_FOC_SNRi.png and b/plots/3C405_x136060/3C405_FOC_SNRi.png differ diff --git a/plots/3C405_x136060/3C405_FOC_SNRp.png b/plots/3C405_x136060/3C405_FOC_SNRp.png index b720859..23f49e5 100755 Binary files a/plots/3C405_x136060/3C405_FOC_SNRp.png and b/plots/3C405_x136060/3C405_FOC_SNRp.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot.png index db4400f..0c2221d 100755 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_IQU.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_IQU.png index 7971e2b..41d72f7 100755 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_IQU.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_IQU.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_P.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_P.png index 449886d..2d2565f 100755 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_P.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_P.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_P_err.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_P_err.png index bd39e1e..5e45a0e 100755 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_P_err.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_P_err.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_P_flux.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_P_flux.png new file mode 100644 index 0000000..ffc7c20 Binary files /dev/null and b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_P_flux.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_SNRi.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_SNRi.png index 4530fb1..954e34c 100755 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_SNRi.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_SNRi.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_SNRp.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_SNRp.png index 77c47e8..051820d 100755 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_SNRp.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM050_rot_SNRp.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot.png index 5677461..d1fb6b1 100755 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_P.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_P.png index 44f4f69..528074c 100755 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_P.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_P.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_P_err.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_P_err.png index 43918b7..048e1d2 100755 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_P_err.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_P_err.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_P_flux.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_P_flux.png index fe450ad..8a446ac 100755 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_P_flux.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_P_flux.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_SNRi.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_SNRi.png index a181389..1079743 100755 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_SNRi.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_SNRi.png differ diff --git a/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_SNRp.png b/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_SNRp.png index a6cd78d..011ffa1 100755 Binary files a/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_SNRp.png and b/plots/3C405_x136060/3C405_FOC_combine_FWHM100_rot_SNRp.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC.png b/plots/CygnusA_x43w0/CygnusA_FOC.png index b56fc9a..257c1a0 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC.png and b/plots/CygnusA_x43w0/CygnusA_FOC.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_IQU.png b/plots/CygnusA_x43w0/CygnusA_FOC_IQU.png index 74acf83..268e14c 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_IQU.png and b/plots/CygnusA_x43w0/CygnusA_FOC_IQU.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_P.png b/plots/CygnusA_x43w0/CygnusA_FOC_P.png index 7477111..b5d4b99 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_P.png and b/plots/CygnusA_x43w0/CygnusA_FOC_P.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_P_err.png b/plots/CygnusA_x43w0/CygnusA_FOC_P_err.png index 9ce6c26..00cce7d 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_P_err.png and b/plots/CygnusA_x43w0/CygnusA_FOC_P_err.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_P_flux.png b/plots/CygnusA_x43w0/CygnusA_FOC_P_flux.png new file mode 100644 index 0000000..c2e9ca2 Binary files /dev/null and b/plots/CygnusA_x43w0/CygnusA_FOC_P_flux.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_SNRi.png b/plots/CygnusA_x43w0/CygnusA_FOC_SNRi.png index 414169c..758724b 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_SNRi.png and b/plots/CygnusA_x43w0/CygnusA_FOC_SNRi.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_SNRp.png b/plots/CygnusA_x43w0/CygnusA_FOC_SNRp.png index 9566fae..6320114 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_SNRp.png and b/plots/CygnusA_x43w0/CygnusA_FOC_SNRp.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot.png index 4af624d..6162c78 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot.png and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_IQU.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_IQU.png index b4fb0e0..cb89369 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_IQU.png and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_IQU.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_P.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_P.png index dcb95be..e860145 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_P.png and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_P.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_P_err.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_P_err.png index 8d39c06..4c269f0 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_P_err.png and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_P_err.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_P_flux.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_P_flux.png new file mode 100644 index 0000000..b4bdf5a Binary files /dev/null and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_P_flux.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_SNRi.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_SNRi.png index 578a9b7..97c1474 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_SNRi.png and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_SNRi.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_SNRp.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_SNRp.png index 33a19b4..5d0d094 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_SNRp.png and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM033_rot_SNRp.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot.png index 3758d20..9eb0b19 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot.png and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_P.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_P.png index f6d0d03..b415e40 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_P.png and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_P.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_P_err.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_P_err.png index 5120b13..20ba1fe 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_P_err.png and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_P_err.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_P_flux.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_P_flux.png index 7d6e920..239c8e6 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_P_flux.png and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_P_flux.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_SNRi.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_SNRi.png index d268160..5c4f213 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_SNRi.png and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_SNRi.png differ diff --git a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_SNRp.png b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_SNRp.png index f94a8cf..f8e1166 100755 Binary files a/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_SNRp.png and b/plots/CygnusA_x43w0/CygnusA_FOC_combine_FWHM066_rot_SNRp.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC.png b/plots/MKN463_x2rp030/MKN463_FOC.png index 44e2ab1..70ef0f0 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC.png and b/plots/MKN463_x2rp030/MKN463_FOC.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_IQU.png b/plots/MKN463_x2rp030/MKN463_FOC_IQU.png index 38f195f..c4b3bf5 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_IQU.png and b/plots/MKN463_x2rp030/MKN463_FOC_IQU.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_P.png b/plots/MKN463_x2rp030/MKN463_FOC_P.png index 7e5471d..a13df09 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_P.png and b/plots/MKN463_x2rp030/MKN463_FOC_P.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_P_err.png b/plots/MKN463_x2rp030/MKN463_FOC_P_err.png index 10608f5..c999e6a 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_P_err.png and b/plots/MKN463_x2rp030/MKN463_FOC_P_err.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_P_flux.png b/plots/MKN463_x2rp030/MKN463_FOC_P_flux.png new file mode 100644 index 0000000..7248a34 Binary files /dev/null and b/plots/MKN463_x2rp030/MKN463_FOC_P_flux.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_SNRi.png b/plots/MKN463_x2rp030/MKN463_FOC_SNRi.png index a27908e..cdc31e9 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_SNRi.png and b/plots/MKN463_x2rp030/MKN463_FOC_SNRi.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_SNRp.png b/plots/MKN463_x2rp030/MKN463_FOC_SNRp.png index 4df8078..08c5099 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_SNRp.png and b/plots/MKN463_x2rp030/MKN463_FOC_SNRp.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot.png index f4c1fb7..7528e87 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot.png and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_IQU.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_IQU.png index 4ee235d..7980562 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_IQU.png and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_IQU.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_P.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_P.png index e1d9baa..a7ecad2 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_P.png and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_P.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_P_err.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_P_err.png index 444c506..86be47c 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_P_err.png and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_P_err.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_P_flux.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_P_flux.png new file mode 100644 index 0000000..a58a5bc Binary files /dev/null and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_P_flux.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_SNRi.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_SNRi.png index 32b23f0..e4bb19c 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_SNRi.png and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_SNRi.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_SNRp.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_SNRp.png index 17ea8f8..8b6e686 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_SNRp.png and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM010_rot_SNRp.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot.png index 858f19c..2445923 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot.png and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_P.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_P.png index 16655ec..e20cddf 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_P.png and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_P.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_P_err.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_P_err.png index 329aea8..364a22a 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_P_err.png and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_P_err.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_P_flux.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_P_flux.png index e8e98cc..5648e36 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_P_flux.png and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_P_flux.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_SNRi.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_SNRi.png index 0c4ad2c..3eba89f 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_SNRi.png and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_SNRi.png differ diff --git a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_SNRp.png b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_SNRp.png index 3654985..2676a1d 100755 Binary files a/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_SNRp.png and b/plots/MKN463_x2rp030/MKN463_FOC_combine_FWHM020_rot_SNRp.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot.png index 801db7b..def3b5f 100755 Binary files a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot.png and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_IQU.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_IQU.png index 05b9da3..a69a9a3 100755 Binary files a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_IQU.png and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_IQU.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_P.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_P.png index 7ffae4a..630b3c7 100755 Binary files a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_P.png and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_P.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_P_err.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_P_err.png index 73df66f..f692725 100755 Binary files a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_P_err.png and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_P_err.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_P_flux.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_P_flux.png new file mode 100644 index 0000000..9e22e77 Binary files /dev/null and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_P_flux.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_SNRi.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_SNRi.png index 915f081..882702a 100755 Binary files a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_SNRi.png and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_SNRi.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_SNRp.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_SNRp.png index dc8f6d3..56aa772 100755 Binary files a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_SNRp.png and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM010_rot_SNRp.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot.png index 1266c8b..ca51df8 100755 Binary files a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot.png and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_P.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_P.png index a2efaf7..3b53cd8 100755 Binary files a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_P.png and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_P.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_P_err.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_P_err.png index ee3ef62..ddd0ea7 100755 Binary files a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_P_err.png and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_P_err.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_P_flux.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_P_flux.png index 12a802a..775f29b 100755 Binary files a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_P_flux.png and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_P_flux.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_SNRi.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_SNRi.png index c469aea..0f8897f 100755 Binary files a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_SNRi.png and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_SNRi.png differ diff --git a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_SNRp.png b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_SNRp.png index 81cb396..7b5b21a 100755 Binary files a/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_SNRp.png and b/plots/NGC1068_x14w010/NGC1068_FOC_combine_FWHM020_rot_SNRp.png differ diff --git a/src/lib/Derivative.cpp b/src/lib/Derivative.cpp deleted file mode 100755 index 10a290c..0000000 --- a/src/lib/Derivative.cpp +++ /dev/null @@ -1,171 +0,0 @@ -// 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; -} - - - -