¿Error in the sampling of positions?

[AsierLambarri]

I was playing around to see how to create a Plummer profile with the limepy class, and performing some checks on some functions i had written. I decided to plot the density profile of the limepy-Plummer i created, sample 1E5 particles from it using the sample class, and then plot the numerically reconstructed profile. I have noticed the following artifact where the density sort of oscillates around the true value given by limepy, as it can be seen in the image bellow:

Sampling 1E5 particles:

I have tried with other numbers of particles and the problem seems to persist:

Sampling 1E6 (top) and 1E4 (bottom) particles:

The oscillating trend can be seen even with the poissonian errorbars. Strangely enough, the problem doesn't seem to exist in the sampling on velocities (i won't attach the plots not to make the issue even longer).

I have tried various binning-rules: scott, friedman-diaconis and sqrt and the problem persists in all of them. In fact, thinner binning and/or larger samples make the issue even more noticeable. With thinner binning one would expect the reconstructed profile to be noisier, but the oscillations still bleed through significantly.

I have found that those oscillations depend on the parameter ode_atol: for smaller absolute tolerances the socillations are smaller, but are still present. Values of ode_atol=1E-15 seem to virtually fix the problem. More fundamentally, the issue may arise from how the positions are sampled (using the inverse CDF method and interpolating r values, if im not mistaken).

Code for reporoducing the error:

from limepy import limepy, sample, spes

k = limepy(
    phi0=0.01,
    M=1E6,
    g=3.49,
    rh=1.3047,
    G=4.300917270038e-06,
    project=True
)

samp = sample(
    k,
    N=10000,
    verbose=True
)

coords = np.vstack((samp.x, samp.y, samp.z)).T
masses = samp.m
vels = np.vstack((samp.vx, samp.vy, samp.vz)).T
radii = np.linalg.norm(coords, axis=1)

bins = 10 ** np.histogram_bin_edges(np.log10(radii), bins="scott")

dens = density_profile(
    unyt_array(coords, 'kpc'),
    unyt_array(masses, 'Msun'),
    bins=bins,
    center=unyt_array([0,0,0], 'kpc')
)

plt.loglog(k.r, k.rho, ls="--")
plt.errorbar(dens['r'], dens['rho'], yerr=dens['e_rho'], ls=":")

plt.ylim(1E-1,)
plt.xlim(5E-3,30)

and the density_profile function:

def density_profile(pos, 
                    mass, 
                    center = None, 
                    bins = None
                   ):
    if center is not None:
        pass
    else:
        center = np.median(pos, axis=0)
        radii = np.linalg.norm(pos - center, axis=1)
        maskcen = radii < 0.5*radii.max()
        center = np.average(pos[maskcen], axis=0, weights=mass[maskcen])
        
        
    coords = pos - center
    radii = np.linalg.norm(coords, axis=1)
    if bins is not None:
        mi, redges, bn = binned_statistic(radii, mass, statistic="sum", bins=bins)
        npbin, _, _ = binned_statistic(radii, mass, statistic="count", bins=bins)

    else:
        mi, redges, bn = binned_statistic(radii, mass, statistic="sum", bins=np.histogram_bin_edges(radii))
        npbin, _, _ = binned_statistic(radii, mass, statistic="count", bins=np.histogram_bin_edges(radii))
    
    redges = redges * coords.units
    mi = mi * mass.units

    if pos.shape[1] == 3:
        volumes = 4 * np.pi / 3 * ( redges[1:]**3 - redges[:-1]**3 )
    if pos.shape[1] == 2:
        volumes = np.pi * ( redges[1:]**2 - redges[:-1]**2 )

    
    rcoords = (redges[1:] + redges[:-1])/2
    dens = mi / volumes
    npbin[npbin == 0] = 1E20
    error = dens / np.sqrt(npbin)
    
    return_dict = {'r': rcoords,
                   'rho': dens,
                   'e_rho': error,
                   'm_enc': mi,
                   'center': center,
                   'dims': pos.shape[1]
                  }
    return return_dict

I am using python3.12, numpy1.26, scipy1.13, unyt3.0.3 and limepy1.1.