Measuring the power spectrum of a generated 3D Gaussian random field (with a specified power spectrum)

I am trying to generate a 3D gaussian random field with a power spectrum of P(k) = 1/k^2, and then measure the power spectrum of the generated field as a consistency check (the measured power spectrum should of course match the analytic one, P(k) = 1/k^2). To generate the 3D field, I used the code https://github.com/cphyc/FyeldGenerator/blob/master/FyeldGenerator/core.py which I found in Creating a 2D Gaussian random field from a given 2D variance. To measure the power spectrum, I used the code in https://github.com/nualamccullagh/zeldovich-bao/blob/master/spatial_stats.py.

Here is the code I used:

  1 import matplotlib.pyplot as plt
  2 import numpy as np
  3 import sys
  4 np.set_printoptions(threshold=sys.maxsize)
  5 import six
  6 import scipy.stats as stats
  7 
  8 #Define global variables
  9 ndim = 3
 10 ngrid = boxsize = 128
 11 n = 2 
 12 A = 1
 13 shape = (ngrid, ngrid, ngrid)
 14 
 15 def generate_field(statistic, power_spectrum, shape, unit_length=1,
 16                                    fft=np.fft, fft_args=dict()):
 17         """                        
 18         Generates a field given a stastitic and a power_spectrum.
 19         """
 20 
 21         fftfreq = np.fft.fftfreq
 22         rfftfreq = np.fft.rfftfreq
 23         
 24         #Compute the k grid     
 25         all_k = [fftfreq(s, d=unit_length) for s in shape[:-1]] + \
 26                         [rfftfreq(shape[-1], d=unit_length)]
 27                         
 28         kgrid = np.meshgrid(*all_k, indexing='ij')
 29         knorm = np.sqrt(np.sum(np.power(kgrid, 2), axis=0))
 30         fourier_shape = knorm.shape
 31         
 32         fftfield = statistic(fourier_shape)
 33         
 34         power_k = np.where(knorm == 0, 0, np.sqrt(power_spectrum(knorm)))
 35         fftfield *= power_k
 36         
 37         return (fft.irfftn(fftfield), fftfield)
 38 
 39 if __name__ == '__main__':
 40         def Pkgen(n):
 41                 def Pk(k):
 42                         return A*np.power(k, -n)
 43                         
 44                 return Pk
 45                 
 46         def distrib(shape):
 47                 # Build a unit-distribution of complex numbers with random phase
 48                 a = np.random.normal(loc=0, scale=1, size=shape)
 49                 b = np.random.normal(loc=0, scale=1, size=shape)
 50                 return a + 1j * b
 51 
 52 # density is the configuration-space density grid (real-valued, shape = ngrid x ngrid x ngrid)
 53 # nkbins is the number of bins to compute the power spectrum in. It computes it in log-spaced bins in the range (2*Pi/L to Pi*ngrid / L)
 54 def getPk(density, nkbins=100):
 55     #make sure the density has mean 0
 56     density=density-np.mean(density)
 57     ngrid=density.shape[0]
 58 
 59     #Fourier transform of density
 60     deltak=np.fft.rfftn(density)
 61     sk=deltak.shape
 62
 63 
 64     #Square the density in Fourier space to get the 3D power, make sure k=0 mode is 0
 65     dk2=(deltak*np.conjugate(deltak)).astype(np.float)
 66     dk2[0,0,0]=0.0
 67 
 68     #set up k-grid
 69     kmin=2*np.pi/boxsize
 70     kny=np.pi*ngrid/boxsize
 71 
 72     a = np.fromfunction(lambda x,y,z:x, sk).astype(np.float)
 73     a[np.where(a > ngrid/2)] -= ngrid
 74     b = np.fromfunction(lambda x,y,z:y, sk).astype(np.float)
 75     b[np.where(b > ngrid/2)] -= ngrid
 76     c = np.fromfunction(lambda x,y,z:z, sk).astype(np.float)
 77     c[np.where(c > ngrid/2)] -= ngrid
 78     kgrid = kmin*np.sqrt(a**2+b**2+c**2).astype(np.float)
 79 
 80     #now we want to compute the power spectrum which involves averaging over shells in k-space 
 81     #define the k-bins we want to compute the power spectrum in
 82     binedges=np.logspace(np.log10(kmin), np.log10(kny),nkbins)
 83     numinbin=np.zeros_like(binedges)
 84     pk=np.zeros_like(binedges)
 85     kmean=np.zeros_like(binedges)
 86     kgridFlatten=kgrid.flatten()
 87     dk2 = dk2.flatten()
 88     index = np.argsort(kgridFlatten)
 89 
 90     kgridFlatten=kgridFlatten[index]
 91     dk2=dk2[index]
 92     c0=0.*c.flatten()+1.
 93     c0[np.where(c.flatten()==0.)]-=0.5
 94     c0=c0[index]
 95     cuts = np.searchsorted(kgridFlatten,binedges)
 96 
 97     for i in np.arange(0, nkbins-1):
 98         if (cuts[i+1]>cuts[i]):
 99             numinbin[i]=np.sum(c0[cuts[i]:cuts[i+1]])
100             pk[i]=np.sum(c0[cuts[i]:cuts[i+1]]*dk2[cuts[i]:cuts[i+1]])
101             kmean[i]=np.sum(c0[cuts[i]:cuts[i+1]]*kgridFlatten[cuts[i]:cuts[i+1]])
102 
103     wn0=np.where(numinbin>0.)
104     pk=pk[wn0]
105     kmean=kmean[wn0]
106     numinbin=numinbin[wn0]
107 
108     pk/=numinbin
109     kmean/=numinbin
110 
111     pk*= boxsize**3/ngrid**6
112 
113     return kmean, pk, kgrid, kmin, a, b, c
114 
115 
116 #call functions
117 densityRealField = np.real(generate_field(distrib, Pkgen(n), shape)[0])
118 km = getPk(densityRealField)[0]
119 PdensityMeasured = getPk(densityRealField)[1]
120 
121 P_analytic = np.zeros(len(km))
122 #Analytic (input) Power Spectrum
123 for i in range(len(km)):
124         P_analytic[i] = A/(km[i]**(n))
125 
126 #plot both analytic and measured power spectrum
127 plt.clf()
128 line1 = plt.plot(km, P_analytic, color = 'cyan', linestyle = 'dashed', label = r'$P(k) \propto 1/k^{2}$')
129 line2 = plt.plot(km, PdensityMeasured, color = 'magenta', label = r'measured $P(k)$')
130 plt.legend()
131 plt.xscale('log')
132 plt.yscale('log')
133 plt.xlabel("$k$")
134 plt.ylabel("$P(k)$")
135 plt.tight_layout()
136 plt.savefig("P_measured_n=2_A=1.png", dpi = 300, bbox_inches = "tight")
137 plt.show()
 

The problem is that the measured power spectrum does not agree with the analytic (i.e. input) power spectrum, as can be seen in this plot. The shape of the measured power spectrum (magenta line) is correct, but it should lie right on top of the analytic one (cyan line). I need the measured power spectrum to ALWAYS match the analytic one (that is, even if I change ngrid values).

Any help is greatly appreciated.