我目前正在分析来自一些流体动力学模拟的一些空间点模式，并且在计算结构因子时遇到了一些困难，$S(\pmb{k})$，从点的位置和径向分布函数，$g(\pmb{r})$. 我一直在关注结构因素维基百科页面以及关于 arxiv 的这篇论文，但我似乎得到了没有意义的结果。

下面是我的类对象，用于计算所述空间点模式的径向分布函数（为清楚起见，我从代码中省略了帮助数据和断言）。该类接受一个 numpy 数组或 numpy 数组列表，找到非零值的位置，计算每个非零值之间的距离，然后计算径向分布。请注意，在这里，如果数组单元格具有非零值，我们将其视为一个点

import numpy as np
import bisect

class pair_correlation_function():

    def __init__(self,data,annulus_width,boundary=None):

        self.data = [data] if type(data) is np.ndarray else data
        self.dr = annulus_width
        self.boundary = boundary

    def get_positions(self):

        positions = list()

        for item in self.data:

            rows,columns = np.nonzero(item)
            positions.append(np.array(list(zip(rows,columns))))

        return positions

    def RDF(self):

        positions = self.get_positions()

        Lx,Ly = np.shape(self.data[0])
        area = Lx*Ly
        radii = list(x for x in range(int(Lx/(2*self.dr))))

        all_radial_distributions = list()

        for item in positions:

            if self.boundary == 'periodic':

                item_new = np.vstack(np.array([np.abs(item[k]-item[k+1:]) for k in range(len(item)-1)]))
                item_new[item_new > Lx/2] = Lx-item_new[item_new > Lx/2]
                norms = [np.sqrt((position**2).sum(axis=0)) for position in item_new]
                norms = sorted(norms)

            else:

                norms = [np.sqrt(((item[k]-item[k+1:])**2).sum(axis=1)) for k in range(len(item)-1)]
                norms = sorted(np.hstack(norms))

            number_particles = len(item)
            item_rdf = list()

            for r in radii:

                i = bisect.bisect_left(norms,r*self.dr)
                j = bisect.bisect_left(norms,(r+1)*self.dr)
                if i != len(norms) and j != len(norms): particle_count = len(norms[i:j])

                normalisation = (2*r+self.dr)*np.pi*self.dr*number_particles**2
                bin_value     = 2*area*particle_count/normalisation

                item_rdf.append(bin_value)

            all_radial_distributions.append(item_rdf)

        radii = np.array(radii)*self.dr

    return radii,all_radial_distributions

下面是一个示例数组（三角形晶格），可以用来玩（我们在这里假设周期性边界条件，因为我的实际数据来自具有周期性边界的模拟）

tri_period = np.array([[0,0,0,0],[0,1,0,1],[0,0,0,0],[1,0,1,0]])
triangular_lattice = np.tile(tri_period,(16,16))
dr = .025

radii, RDF = pair_correlation_function(triangular_lattice,dr,'periodic').RDF()

和绘图给出

在这里，我们看到了一系列预期的布拉格峰（尽管我不确定为什么峰会分割成平滑、单调的递减带？但我的实际数据并没有这样做并且看起来是正确的。）现在，如果我尝试计算给定位置的结构因子（参见上面的 arxiv 论文，方程$(2)$)

def sf_from_positions(positions,box_dimension):

    sf = list()
    modes = list(x for x in range(1,int(box_dimension)))

    dk = 2*np.pi/box_dimension

    for h in modes:

        k_vec = np.array([1,1])*dk*h
        summation = 0

        for position in positions:

            summation += np.exp(-1j*k_vec.dot(position))

        sf.append(abs(summation)**2/len(positions))

    return sf, modes

rows,columns = np.nonzero(triangular_lattice)
positions = list(zip(rows,columns))
dimension = np.shape(triangular_lattice)[0]
factor, wave_modes = sf_from_positions(positions,dimension)

和绘图给出

看起来不对，因为布拉格峰的傅里叶变换也应该是布拉格峰（请参阅此 arxiv 论文，第 12 页）。

所以，我的第一个问题是，我在计算结构因子时出了什么问题？谁能看到我做错了什么？其次，我想知道如何通过公式整合径向分布函数来获得结构因子

$$S(\pmb{k}) = 1 + 4 \pi \rho \int_{\mathbb{R^{n}}} [g(\pmb{r})-1] e^{-i \pmb{k} \cdot \pmb{r}} d \pmb{r}$$ 我尝试在 scipy 中使用 cumtrapz 模块，但它给了我一个负数（不可能）并且随着波包增加而振荡的结构因子。

对不起，如果这个问题过于牵涉。谢谢你的帮助。