Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
from pylab import *
import astropy.io.fits as fits
import xqml
import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
nside = 8
lmax = 3*nside-1
Slmax = 3*nside-1
deltal = 1
clth = read_cl( 'planck_base_planck_2015_TTlowP.fits')
lth = arange(2,lmax+1)
ellbins = arange(2,lmax+2,deltal)
ellbins[-1] = lmax+1
P, Q, ell, ellval = xqml.GetBinningMatrix(ellbins, lmax)
nbins=len(ellbins)-1
#create mask
t,p=hp.pix2ang( nside, range(hp.nside2npix(nside)))
mask = np.ones(hp.nside2npix(nside),bool)
mask[abs(90-rad2deg(t)) < 10] = False
ip=arange(hp.nside2npix(nside))
ipok=ip[mask]
npix=len(ipok)
fwhm = 0.5
bl=gauss_beam(deg2rad(fwhm), lmax=Slmax+1)
allStoke, der, ind = xqml.getstokes(polar=True,temp=False,EBTB=False)
nder=len(der)
muKarcmin = 0.1
###################################### Compute ds_dcb ######################################
Pl, S = xqml.compute_ds_dcb(ellbins,nside,ipok,bl,clth,Slmax,polar=True,temp=False, EBTB=False, pixwining=True, timing=True, MC=False)
Pl = Pl.reshape((nder)*(np.shape(Pl)[1]), 2*npix, 2*npix)
###################################### Construct El, Wb ######################################
pixvar = 2*xqml.muKarcmin2var(muKarcmin, nside)
varmap = ones((2*npix))*pixvar
NoiseVar = np.diag(varmap)
C = S + NoiseVar
invC = inv(C)
#invS = inv(S)
El = xqml.El(invC, invC, Pl)
Wb = xqml.CrossWindowFunction(El, Pl)
invW = inv(Wb)
###################################### Compute spectra ######################################
cmb = hp.synfast( clth, nside, fwhm=deg2rad(fwhm), pixwin=True, new=True, verbose=False)
noise = (randn(len(varmap)) * varmap**.5).reshape(2,-1)
dm = cmb[1:,mask] + noise
yl = xqml.yQuadEstimator(dm.ravel(), dm.ravel(), El)
cl = xqml.ClQuadEstimator(invW, yl).reshape(nder,-1)
plot( clth.transpose()[:lmax+1,:])
plot( ellval, cl.transpose())
G = xqml.CrossGisherMatrix(El, C)
V = xqml.CovAB(invW, G)
###################################### Construct MC ######################################
allcl = []
for n in range(100):
progress_bar( n, 100)
cmb = hp.synfast( clth, nside, fwhm=deg2rad(fwhm), pixwin=True, new=True, verbose=False)
dm = cmb[1:,mask] + (randn(2*npix)*sqrt(varmap)).reshape( 2, npix)
yl = xqml.yQuadEstimator(dm.ravel(), dm.ravel(), El)
allcl.append( xqml.ClQuadEstimator(invW, yl).reshape(nder,-1))
fits.writeto( "allcl_Slmax%d.fits" % Slmax, array(allcl))
figure()
subplot( 2,1,1)
plot( lth, clth.transpose()[lth,1:3], 'k')
plot( ellval, mean( allcl,0).transpose(), 'r')
plot( ellval, mean( allcl,0).transpose() + std( allcl,0).transpose(), 'r--')
plot( ellval, mean( allcl,0).transpose() - std( allcl,0).transpose(), 'r--')
subplot( 2,1,2)
cosmic = sqrt( 2./(2*lth+1))/mean(mask) * clth[1:3,lth]
plot( lth, cosmic.transpose(), 'k')
plot( ellval, std( allcl,0).transpose(), 'r')
plot( ellval, sqrt(diag(V)).reshape(nder,-1).transpose(), 'b')