Tests for code that computes two-point correlation function of galaxies

Question

I'm currently writing a code to compute the two-point correlation function of a set of galaxies. As I work with a loooot of galaxies, it was suggested to me not to do it in real space, but to compute the power spectrum in Fourier space and get the correlation from the power spectrum, as is explained e.g. in appendix B in this paper (arXiv:1507.01948).

What I need are some very simple test cases with known results to test my code at each step in order to ensure that I really get what I want. Can anybody point me to any such simple tests? Or at least to a already existing code (preferably in Python or Fortran) that does exactly this, so I can cross-compare the results?

Answer 1

A very simple example is to fill a 3d density with a regular plane wave, which should result in a single peak in the power spectrum. Here is a simple python scripts that does exactly that.

Computing the two point correlation function should be only one inverse fourier transform of the Power spectrum away. (Depending on which FFT library you use, you might need to renormalize the results. FFTW and numpy.fft for example have unnormalized fourier transforms: $F^{-1}[F[f(x)] = Nf(x) $, where $N$ is the number of samples.)

#!/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()