Pruebas de código que calcula la función de correlación de dos puntos de las galaxias

Actualmente estoy escribiendo un código para calcular la función de correlación de dos puntos de un conjunto de galaxias. Como trabajo con un montón de galaxias, me sugirieron que no lo hiciera en el espacio real, sino que calculara el espectro de potencia en el espacio de Fourier y obtuviera la correlación del espectro de potencia, como se explica, por ejemplo, en el apéndice B de este artículo . (arXiv:1507.01948) .

Lo que necesito son algunos casos de prueba muy simples con resultados conocidos para probar mi código en cada paso para asegurarme de que realmente obtengo lo que quiero. ¿Alguien puede señalarme alguna prueba tan simple? ¿O al menos a un código ya existente (preferiblemente en Python o Fortran) que hace exactamente esto, para que pueda comparar los resultados?

Para probar su código, presumiblemente puede usar cualquier par de espectros de potencia, en lugar de usar alguna fuente de datos astrofísicos. Intente expandir su búsqueda a "introducción a la correlación espectral de potencia" o algo así.
Lo intenté, pero no puedo encontrar ninguno para campos (estoy trabajando con campos de densidad 3D), los ejemplos suelen ser todos 1D...

Respuestas (1)

Un ejemplo muy simple es llenar una densidad 3D con una onda plana regular, lo que debería dar como resultado un solo pico en el espectro de potencia. Aquí hay un script de Python simple que hace exactamente eso.

El cálculo de la función de correlación de dos puntos debe estar a solo una transformada de Fourier inversa del espectro de potencia. (Dependiendo de la biblioteca FFT que utilice, es posible que deba volver a normalizar los resultados. FFTW y numpy.fft, por ejemplo, tienen transformadas de Fourier no normalizadas: F − 1 [ F [ F ( X ) ] = norte F ( X ) , donde norte es el número de muestras).

#!/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()
Pruebas de código que calcula la función de correlación de dos puntos de las galaxias
RespuestasAquíPolítica de privacidad1y, 0v