我正在使用HEALPy的anafast从普朗克地图中提取Cl,或者是数据地图
1
,或者模拟地图。似乎即使在我使用了它们的强度普通面膜之后
2
,我得到的功率谱大约是他们释放的五倍
3
..。
# extract Cl's
filename = 'path/to/COM_CMB_IQU-commander_2048_R3.00_full.fits'
test_map = read_map(filename)
path = 'path/to/COM_Mask_CMB-common-Mask-Int_2048_R3.00.fits'
mask = hp.read_map(path)
map_masked = hp.ma(test_map)
map_masked.mask = np.logical_not(mask)
test_cls_meas_frommap = anafast(map_masked, lmax=3000)
# load up Planck meas powerspectrum
path = '/path/to/powerspectrum'
T2 = 10**12*2.7255**2 # in muK^2
datalist = {
'tt': 'COM_PowerSpect_CMB-TT-binned_R3.01.txt',
'ee': 'COM_PowerSpect_CMB-EE-binned_R3.02.txt',
'te': 'COM_PowerSpect_CMB-TE-binned_R3.02.txt'
}
targ = os.path.join(path, datalist['tt'])
res = cmb.load_meas(targ)
ll_meas = res[:, 0]
test_cls_meas = res[:, 1]/ll_meas/(ll_meas+1)*2.*np.pi/T2
# output
plt.subplots()
plt.plot(ll_meas, ll_meas*(ll_meas+1.)*test_cls_meas*T2/2./np.pi, '--', alpha=0.6, label='Planck 2018 PS release')
plt.plot(ll, ll*(ll+1.)*test_cls_meas_frommap*T2/2./np.pi, '--', alpha=0.6, label='Planck 2018 PS from Data Map')
plt.xlabel(r'$\ell$')
plt.ylabel(r'$D_\ell$')
#plt.xscale('log')
plt.legend(loc='best')另一方面,如果我自己使用synfast合成地图,然后使用anafast提取功率谱,我可以验证我是否获得了输入功率谱。我想知道与普朗克方法相比,是否有任何潜在的陷阱可能导致功率谱计算的不匹配?
数据来源:
1
数据映射:(wget -O COM
_
CMB
_
IQU-指挥官
_
2048
_
R3.00
_
完整"pla.esac.esa.int/pla-sl/data-action?MAP.MAP
_
OID=13470")
2
掩码映射:(wget -O COM
_
遮罩
_
CMB-common-Mask-Int
_
2048
_
R3.00.fits“
http://pla.esac.esa.int/pla/aio/product-action?MAP.MAP
_
ID=COM
_
遮罩
_
CMB-common-Mask-Int
_
2048
_
R3.00.fits“
)
3
官方功率谱:(wget -O COM
_
PowerSpect
_
CMB-TT-binned
_
R3.01.txt“
http://pla.esac.esa.int/pla/aio/product-action?COSMOLOGY.FILE
_
ID=COM
_
PowerSpect
_
CMB-TT-binned
_
R3.01.txt“
)
发布于 2021-03-01 02:10:15
您的计算中有几个问题:
入库的功率谱已经是D
_
哦,你不应该把它乘以ell(ell+1)因子,也不应该乘以T2
当你在切割天空上计算光谱时,你必须除以天空分数。
所以这应该是可行的:
cmb_binned_spectrum = np.loadtxt('COM_PowerSpect_CMB-TT-binned_R3.01.txt')
k2muK = 1e6
plt.plot(cmb_binned_spectrum[:,0], cmb_binned_spectrum[:,1], '--', alpha=1, label='Planck 2018 PS release')
plt.plot(ll, ll*(ll+1.)*test_cls_meas_frommap*k2muK**2/2./np.pi / sky_fraction, '--', alpha=0.6, label='Planck 2018 PS from Data Map')
plt.xlabel(r'$\ell$')
plt.ylabel(r'$D_\ell~[\mu K^2]$')
plt.grid()
plt.legend(loc='best')我解释了它,并完成了这个笔记本中的所有步骤:
https://zonca.dev/2021/02/compute-planck-spectra-healpy-anafast.html
https://stackoverflow.com/questions/66407984
复制相似问题