Update compilation and script test

examples/script_xpol.py

@@ -4,7 +4,7 @@ import healpy as hp
 from numpy import *
 import matplotlib.pylab as plt
 
-nside=128
+nside=256
 lmax=3*nside-1
 tagnames = ['TT','EE','BB','TE','TB','EB']
 
@@ -30,6 +30,7 @@ xp = xpol.Xpol( mask, binning, verbose=True)
 #compute spectra from map dT
 pcl,cl = xp.get_spectra( dT)
 
+
 #plot result
 plt.figure()
 plt.subplot( 2,1,1)
@@ -43,39 +44,42 @@ plt.axhline( 1, color='k')
 
 
 #small MonteCarlo
+nsimu = 100
 xp = xpol.Xpol( mask, binning)
 allcl = []
-for i in range(100):
+for i in range(nsimu):
     print(i)
     dT = hp.synfast( clth, nside, new=True, lmax=3*nside-1, pixwin=True, verbose=False)
     pcl,cl = xp.get_spectra( dT, pixwin=True)
     allcl.append(cl)
 
 
+#figure
 cb = binning.bin_spectra(clth)
 cb = list(cb)+list([cb[0]*0.,cb[0]*0.])
-## cosmic = sqrt( 2./(2*l+1))/mean(mask) * cb
-## cosmic[3] = sqrt( 1./(2*l+1)*(cb[0]*cb[1]+cb[3]*cb[3]))/xpol.fsky(mask) 
 cosmic = xpol.sample_variance( lb, cb, fsky=xpol.fsky(mask))
 
 tag=1
 plt.figure()
 for tag in range(6):
     plt.subplot( 3,6,0*6+tag+1)
-    plt.plot( lb, lb*(lb+1)/2/pi*cb[tag], 'k')
-#    plt.semilogx()
-    plt.plot( lb, lb*(lb+1)/2/pi*mean(allcl,0)[tag], 'r')
-    plt.plot( lb, lb*(lb+1)/2/pi*(mean(allcl,0)-std(allcl,0))[tag], 'r--')
-    plt.plot( lb, lb*(lb+1)/2/pi*(mean(allcl,0)+std(allcl,0))[tag], 'r--')
+    plt.plot( lb, lb*(lb+1)/2/pi*cb[tag]*1e12, 'k')
+    plt.plot( lb, lb*(lb+1)/2/pi*mean(allcl,0)[tag]*1e12, 'r')
+    plt.plot( lb, lb*(lb+1)/2/pi*(mean(allcl,0)-std(allcl,0))[tag]*1e12, 'r--')
+    plt.plot( lb, lb*(lb+1)/2/pi*(mean(allcl,0)+std(allcl,0))[tag]*1e12, 'r--')
     plt.title( tagnames[tag])
-    plt.ylabel( "$\ell(\ell+1)/2\pi\,C_\ell$")
+    if tag == 0: plt.ylabel( "$\ell(\ell+1)/2\pi\,C_\ell\,[\mu K^2]$")
     plt.subplot( 3,6,1*6+tag+1)
-    plt.plot( lb, mean(allcl,0)[tag]/cb[tag], 'r')
-    plt.axhline( 1, color='k')
+    plt.plot( lb, (mean(allcl,0)[tag]-cb[tag])/(std(allcl,0)[tag]/sqrt(nsimu)), 'r')
+    plt.ylim( -4,4)
+    plt.axhline( 0, color='k')
+    if tag == 0: plt.ylabel( "$R_\ell [\sigma]$")
     plt.subplot( 3,6,2*6+tag+1)
-    plt.plot( lb, lb*(lb+1)/2/pi*std(allcl,0)[tag], 'r')
-    plt.plot( lb, lb*(lb+1)/2/pi*cosmic[tag], 'k')
-    plt.ylabel( "$\ell(\ell+1)/2\pi\,\sigma_\ell$")
+    plt.plot( lb, lb*(lb+1)/2/pi*std(allcl,0)[tag]*1e12, 'r')
+    plt.plot( lb, lb*(lb+1)/2/pi*cosmic[tag]*1e12, 'k')
+    if tag == 0: plt.ylabel( "$\ell(\ell+1)/2\pi\,\sigma_\ell\,[\mu K^2]$")
+
+
 
 
 

makefile_gnu

@@ -4,12 +4,12 @@ HEALPIX=/sps/planck/Users/tristram/Soft/Xpol/software/libs/Healpix_3.70
 CFITSIO=/sps/planck/Users/tristram/Soft/Xpol/software/libs/
 S2HAT=/sps/planck/Users/tristram/Soft/Xpol/software/libs/s2hat
 
-export INC_DIR = -I${S2HAT_DIR} -I${HEALPIX_DIR}/include -I${CFITSIO_DIR}/include
+export INC_DIR = -I${MKLROOT}/include -I${S2HAT_DIR} -I${HEALPIX_DIR}/include -I${CFITSIO_DIR}/include
 export LIB_DIR = \
           -L${S2HAT_DIR} -ls2hat_hlpx_g \
           -Wl,-R${HEALPIX_DIR}/lib -L${HEALPIX_DIR}/lib -lchealpix -lhealpix -lsharp -lm \
           -Wl,-R${CFITSIO_DIR} -L${CFITSIO_DIR} -lcfitsio -lm \
-	  -Wl,-R${MKLROOT} -L${MKLROOT} -lmkl_blacs_openmpi_lp64 -lmkl_gf_lp64 -lmkl_core -lmkl_gnu_thread -lmkl_def -liomp5 -ldl -lpthread -lm
+	  -Wl,-R${MKLROOT}/lib/intel64 -L${MKLROOT}/lib/intel64 -lmkl_blacs_openmpi_lp64 -lmkl_gf_lp64 -lmkl_core -lmkl_gnu_thread -lmkl_def -ldl -lpthread -lm
 
 export LIB = 
 
@@ -19,7 +19,8 @@ export F77 = mpif77
 
 export CFLAG = -O3 -DMKL -DHEALPIXDATA=\"${HEALPIX_DIR}/data\" -Wall -Dx86_64 -DOPEN_MPI
 export FFLAG = -O3 -fno-second-underscore
-export LFLAG = -O3 -lmpi -lmpi_mpifh -lgfortran -lgomp -fopenmp #-lifcore -lintlc -lifport -limf -lirc -ldl -lsvml
+export LFLAG = -O3 -lmpi -lmpi_mpifh -lgfortran -lgomp -fopenmp
+#-liomp5 -lifcore -lintlc -lifport -limf -lirc -ldl -lsvml
 ################################################################################################