Тесты для кода, вычисляющего двухточечную корреляционную функцию галактик

В настоящее время я пишу код для вычисления функции двухточечной корреляции набора галактик. Поскольку я работаю с большим количеством галактик, мне было предложено не делать это в реальном пространстве, а вычислять спектр мощности в пространстве Фурье и получать корреляцию из спектра мощности, как объяснено, например, в приложении B к этой статье . (архив: 1507.01948) .

Что мне нужно, так это несколько очень простых тестов с известными результатами для проверки моего кода на каждом этапе, чтобы убедиться, что я действительно получаю то, что хочу. Может ли кто-нибудь указать мне на какие-либо такие простые тесты? Или, по крайней мере, к уже существующему коду (желательно на Python или Fortran), который делает именно это, чтобы я мог сравнить результаты?

Предположительно, для проверки вашего кода вы можете использовать любые пары спектров мощности, а не использовать какой-либо источник астрофизических данных. Попробуйте расширить поиск до «введение в спектральную корреляцию мощности» или что-то в этом роде.
Я пытался, но не могу найти поля для полей (я работаю с полями 3D-плотности), примеры обычно все 1D...

Ответы (1)

Очень простой пример — заполнить трехмерную плотность регулярной плоской волной, что должно привести к одному пику в спектре мощности. Вот простой скрипт Python, который делает именно это.

Вычисление двухточечной корреляционной функции должно быть на расстоянии всего одного обратного преобразования Фурье спектра мощности. (В зависимости от того, какую библиотеку БПФ вы используете, вам может потребоваться перенормировать результаты. Например, FFTW и numpy.fft имеют ненормализованные преобразования Фурье: Ф − 1 [ Ф [ ф ( Икс ) ] знак равно Н ф ( Икс ) , куда Н это количество образцов.)

#!/usr/bin/python3

import numpy as np
import matplotlib.pyplot as plt


#==================================
def main():
#==================================


    nc = 128                # define how many cells your box has
    boxlen = 50.0           # define length of box
    Lambda = boxlen/4.0     # define an arbitrary wave length of a plane wave
    dx = boxlen/nc          # get size of a cell

    # create plane wave density field
    density_field = np.zeros((nc, nc, nc), dtype='float')
    for x in range(density_field.shape[0]):
        density_field[x,:,:] = np.cos(2*np.pi*x*dx/Lambda)

    # get overdensity field
    delta = density_field/np.mean(density_field) - 1

    # get P(k) field: explot fft of data that is only real, not complex
    delta_k = np.abs(np.fft.rfftn(delta).round())
    Pk_field =  delta_k**2

    # get 3d array of index integer distances to k = (0, 0, 0)
    dist = np.minimum(np.arange(nc), np.arange(nc,0,-1))
    dist_z = np.arange(nc//2+1)
    dist *= dist
    dist_z *= dist_z
    dist_3d = np.sqrt(dist[:, None, None] + dist[:, None] + dist_z)

    # get unique distances and index which any distance stored in dist_3d 
    # will have in "distances" array
    distances, _ = np.unique(dist_3d, return_inverse=True)

    # average P(kx, ky, kz) to P(|k|)
    Pk = np.bincount(_, weights=Pk_field.ravel())/np.bincount(_)

    # compute "phyical" values of k
    dk = 2*np.pi/boxlen
    k = distances*dk

    # plot results
    fig = plt.figure(figsize=(9,6))
    ax1 = fig.add_subplot(111)
    ax1.plot(k, Pk, label=r'$P(\mathbf{k})$')

    # plot expected peak:
    # k_peak = 2*pi/lambda, where we chose lambda for our planar wave earlier
    ax1.plot([2*np.pi/Lambda]*2, [Pk.min()-1, Pk.max()+1], label='expected peak')
    ax1.legend()
    plt.show()








#==================================
if __name__ == "__main__":
#==================================

    main()
Тесты для кода, вычисляющего двухточечную корреляционную функцию галактик
СпросиСетьПолитика конфиденциальности1y, 0v