我的三维离散功率谱计算代码有什么问题?(3D power spectrum)

为了计算一个256^3的密度分布的功率谱,我写了下面的代码,但计算得到的图像与应该出现的图像不一样,我想问一下是哪里出了问题

import numpy as np
from numpy import fft
import scipy.stats as stats
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
import pandas as pd
import math


#打开密度分布文件
with open('D:\\science data\\1\\dm_dens_resamp_0_01.dat')as ff:
    dens_field=np.fromfile(ff,
                            dtype='float32').reshape((256,256,256),order='F')
Lbox=75
Ng=256

#Taking the Fourier transform
A=dens_field-1#A是网格化后的相对密度
npix=A.shape[0]
fouier_dens_field=np.fft.fftn(A)
B=fouier_dens_field#B是傅里叶变换后得到的复数组
fourier_amplitudes=np.abs(B)**2
C=fourier_amplitudes#C是B的模方
D=C*(pow(Lbox,3)/pow(Ng,6))

#Constructing a wave vector array
kfreq=np.fft.fftfreq(npix)*npix#k空间波矢分布
kfreq3D=np.meshgrid(kfreq,kfreq,kfreq)#转换成对应于三位密度场数据的三维数组
knrm=np.sqrt(kfreq3D[0]**2+kfreq3D[1]**2+kfreq3D[2]**2)#我们考虑的是它的模
    #由于要将波矢模与傅里叶结果作为两个维度输出图像,
    #因此将上述数组展平
knrm=knrm.flatten()
fourier_amplitudes=D.flatten()

#Creating the power spectrum
kbins=np.arange(0.5,npix//2+1,1.)#用于在k空间中储存振幅
kvals=0.5*(kbins[1:]+kbins[:-1])#放入数据
    #以下计算每个kbin中平均傅里叶振幅(平方)
Abins, _, _ = stats.binned_statistic(knrm,fourier_amplitudes,
                                     statistic = "mean",
                                     bins = kbins)#总方差
Abins *=4.*np.pi/3.*(kbins[1:]**3-kbins[:-1]**3)#还需要乘每个bin的体积

#画图
fig=plt.figure(num=1, figsize=(40, 40),dpi=50)
ax1=fig.add_subplot(2,1,1)
ax1.set_title('256^3 Power Spectral Density')
ax1.set_xlabel('K(x^-1)')
ax1.set_ylabel('y-PSD')
plt.loglog(kvals,Abins)
plt.show()

你好，我是有问必答小助手，非常抱歉，本次您提出的有问必答问题，技术专家团超时未为您做出解答

本次提问扣除的有问必答次数，将会以问答VIP体验卡（1次有问必答机会、商城购买实体图书享受95折优惠）的形式为您补发到账户。

​​​​因为有问必答VIP体验卡有效期仅有1天，您在需要使用的时候【私信】联系我，我会为您补发。