Testet auf Code, der die Zweipunkt-Korrelationsfunktion von Galaxien berechnet

Ein sehr einfaches Beispiel ist das Füllen einer 3D-Dichte mit einer regelmäßigen ebenen Welle, was zu einer einzelnen Spitze im Leistungsspektrum führen sollte. Hier ist ein einfaches Python-Skript, das genau das tut.

Die Berechnung der Zweipunkt-Korrelationsfunktion sollte nur eine inverse Fourier-Transformation des Leistungsspektrums entfernt sein. (Je nachdem, welche FFT-Bibliothek Sie verwenden, müssen Sie möglicherweise die Ergebnisse neu normalisieren. FFTW und numpy.fft zum Beispiel haben nicht normalisierte Fourier-Transformationen:$F^{-1}[F[f(x)] = Nf(x)$, wo$N$ist die Anzahl der Proben.)

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