Computing the structure factor from positions and radial distribution function [closed]

Question

I'm currently analysing some spatial point patterns that come from some fluid dynamics simulations and I'm having some difficulty computing the structure factor, $S(\pmb{k})$, from both the positions of the points and the radial distribution function, $g(\pmb{r})$. I have been following the structure factor Wikipedia page and also this paper on arxiv but I seem to be getting results that don't make sense.

Below is my class object to compute the radial distribution function for said spatial point patterns (I've omitted the help data and assertions from the code for clarity). The class takes in a numpy array or list of numpy arrays, finds the locations of the nonzero values, computes the distances between each nonzero value and then computes the radial distribution. Note that here, we consider an array cell a point if it has a nonzero value

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

Below is a sample array (a triangular lattice) with which to play around (we are assuming periodic boundary conditions here because my actual data comes from a simulation with periodic boundaries)

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()

and plotting gives

Here, we see a series of Bragg peaks as expected (though I'm not sure why the peaks segment into smooth, monotone decreasing bands? My actual data doesn't do this though and looks correct.) Now, if I try to compute the structure factor given the positions (see arxiv paper above, equation $(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)

and plotting gives

which looks off, because the Fourier transform of Bragg peaks should also be Bragg peaks (see this arxiv paper, page 12).

So, my first question is, what is going wrong in my computation of the structure factor? Can anyone see what I have done wrong? And secondly, I was wondering how I would integrate the radial distribution function to get the structure factor, through the formula $$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} $$ I tried to use the cumtrapz module in scipy but it gave me a structure factor that was negative (impossible) and that oscillated with an increasing wave packet.

Sorry if the question is overly involved. Thanks for your help.