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

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?

• To test your code, presumably you can use any pairs of power spectra, rather than using some astrophysical data source. Try expanding your search to "intro to power spectral correlation" or some such. – Carl Witthoft May 29 '18 at 14:14
• I did try, but I can't find any for fields (I'm working with 3d density fields), the examples are usually all 1D... – mivkov May 31 '18 at 9:34

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.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()