У меня есть сетка данных, где строки представляют тета (0, пи), а столбцы представляют фи (0, 2 * пи), и где f (тета, фи) - плотность темной материи в этом месте. Я хотел рассчитать спектр мощности для этого и решил использовать healpy.

Что я не могу понять, так это как отформатировать мои данные, чтобы использовать их. Если бы кто-то мог предоставить код (на python по понятным причинам) или указать мне учебник, это было бы здорово! Я попытался сделать это с помощью следующего кода:

#grid dimensions are Nrows*Ncols (subject to change)
theta = np.linspace(0, np.pi, num=grid.shape[0])[:, None]
phi = np.linspace(0, 2*np.pi, num=grid.shape[1])
nside = 512
print "Pixel area: %.2f square degrees" % hp.nside2pixarea(nside, degrees=True)
pix = hp.ang2pix(nside, theta, phi)
healpix_map = np.zeros(hp.nside2npix(nside), dtype=np.double)
healpix_map[pix] = grid 

Но, когда я пытаюсь выполнить код, чтобы сделать спектр мощности. В частности,

cl = hp.anafast(healpix_map[pix], lmax=1024)

Я получаю эту ошибку:

Ошибка типа: неверное количество пикселей

Если бы кто-нибудь мог указать мне хороший учебник или помочь отредактировать мой код, это было бы здорово.

Дополнительные характеристики: мои данные находятся в массиве 2d np, и я могу изменить numRows / numCols, если потребуется.

Редактировать:

Я решил эту проблему, сначала изменив аргументы anafast на healpix_map. Я также улучшил интервал, сделав свои Nrows * Ncols = 12 * nside * nside. Но мой энергетический спектр все еще дает ошибки. Если у кого-то есть ссылки на хорошую документацию / учебное пособие о том, как рассчитать спектр мощности (состояние тэта / фи арг), это было бы невероятно полезно.