The following code draws a random distribution of particles in 3D. The Manipulate box allows to change a few parameters (number of particles, size of clusters, ...). I would like to define the local number of particles per unit volume $V = \ell^3$, on a variable scale $\ell$, and plot it as level curves (i.e surfaces) in 3D. How can this be done in Mathematica?

galaxies[p_, q_, r_] := Module[
    {pts = RandomReal[{-1, 1}, {p, 3}]},
    pts = Nest[Join[#, {RandomChoice[#] + r RandomReal[{-1, 1}, 3]}] &, pts, q];

graph[p_, q_, r_] := Graphics3D[{RGBColor[{0.5, 0.4, 1.0}], PointSize[0.004], Point[galaxies[p, q, r]]}]

view[p_, q_, r_] := Show[{graph[p, q, r]},
    PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}},
    Boxed -> True,
    Background -> Black,
    ImageSize -> {700, 700},
    SphericalRegion -> True,
    Method -> {"RotationControl" -> "Globe"}

    view[p, q, r],
    {{p, 500, "Number of clusters"}, 1, 1000, 1},
    {{q, 10000, "Number of galaxies"}, 0, 10000, 1},
    {{r, 0.05, "Scale"}, 0, 0.2, 0.001},
    ControlPlacement -> Bottom

The scale $\ell$ should be a variable, from a small value up to the full box size, say : $0.01 \le \ell \le 2$, that could be changed in another Manipulate box. The idea is to plot the density of the particles in the box as a scalar function: $$\tag{1} \rho(x, y, z, \ell) = \frac{N(x, y, z, \ell)}{\ell^3}, $$ where $N(x, y, z, \ell)$ gives the total number of particles inside the cube of volume $V = \ell^3$, as a function of the cube center $\{x, y, z \}$ (with -1 < x, y, z < 1). I really don't know how to start to implement the idea, especially since the sampling volume should always stay inside the whole box of size $\ell_{max} = 2$ (with center located at $\{0, 0, 0 \}$). The sampling near the edges of the whole box may also be a problem (maybe the random distribution should be defined on a box two times larger than the display box).

EDIT: The exact "microscopic" density of the $N$ particles can be defined using Dirac's deltas: $$\tag{2} \rho(x, y, z) = \sum_{k \,=\, 1}^N \delta (x - x_k) \, \delta (y - y_k) \, \delta (z - z_k). $$ This density is 0 everywhere, except at the location of a particle. We could then introduce the averaged density on a volume $\ell^3$ (around a point of coordinates $\{x, y, z\}$), by integration: $$\tag{3} \bar{\rho}_{\ell}(x, y, z) = \frac{1}{\ell^3} \int_{x - \ell/2}^{x + \ell/2} \int_{y - \ell/2}^{y + \ell/2} \int_{z - \ell/2}^{z + \ell/2} \rho(x, y, z) \, dx \, dy \, dz. $$ That density could be called "fine grained" when $\ell$ is small, and "coarse grained" when $\ell$ is close to the size of the whole box. In all cases: $0.0001 \le \ell \le 2$ (since the random distribution above is enclosed in a box of size $2$).

I should have stated this already from the start in my question, sorry about that. The parameter $\ell$ should be a variable in a Manipulate box. The function (3) is what I would like to calculate for my random distribution (code above), and plot as a 3D graphics (level curves or surfaces?). Of course, we should get $\bar{\rho}_1(x, y, z) = N/8$ when $\ell = 2$.

Ideally, the density $\bar{\rho}_{\ell}(x, y, z)$ should be smoothed/interpolated to give some smooth level curves/surfaces in a plot, instead of ugly discontinuous jumps. So how should I do this, with Mathematica?

  • $\begingroup$ Unless I'm missing something, this sounds like a job for SmoothKernelDistribution $\endgroup$
    – Lukas Lang
    Commented Aug 5, 2022 at 16:43
  • $\begingroup$ @LukasLang, hmm, this sounds interesting. Unfortunately, this function isn't implemented in my old Mma 7.0. $\endgroup$
    – Cham
    Commented Aug 5, 2022 at 16:45
  • $\begingroup$ You can probably implement something similar by hand. All the relevant equations are given in the documentation I think. Two options that come to mind are to define a grid and add the contributions of each point into that grid (i.e. to effectively compute the distribution on a predefined grid), or to define a function as the sum of the different contributions (effectively the first equation from the documentation), and let the plotting function do the samplong itself. $\endgroup$
    – Lukas Lang
    Commented Aug 5, 2022 at 16:54
  • $\begingroup$ I might certainly be not understanding the question but why would the average density change dependent on the size of the sampling unit? (Whether that sampling unit be a box or a sphere or some other shape). However, the density would certainly vary more among samples with smaller sized sampling units. $\endgroup$
    – JimB
    Commented Aug 5, 2022 at 17:46
  • 1
    $\begingroup$ @HansOlo I think you bring up a good question for Mathematica Meta. I'm pretty judgmental about the statistical aspects of questions on this site but helping someone with a research project "for free" doesn't bother me so much. $\endgroup$
    – JimB
    Commented Aug 5, 2022 at 19:50

2 Answers 2


The following is an example based on a uniform distribution of points/galaxies, which can be adapted to any desired distribution. Also, some notes:

  1. The code is given here only for educational purposes and several optimizations can be made.
  2. Using Nearest[] is not the most efficient approach to calculate "counts in cells", ie the number of particles in growing spheres of radius $r$ around each galaxy. It is in fact rather inefficient.

First, let's make two populations of points

data = RandomReal[{-2, 2}, {2000, 3}];
datarmin1 = Select[data, #.# <= 1 &];

For this approach (and all similar "counts in cells" methods), one needs either more overall galaxies (given by data) than then ones under study (given by datarmin1) or periodic boundary conditions (eg coming from an N-body sim).

Then we use Nearest[] to count the number $N(r)$ of points around each of the galaxies under study, within the larger sample:

r0 = 0.01;
rmax = 1.01;
dr = 0.05;
j`j = 0; SetSharedVariable[j`j];
ProgressIndicator[Dynamic[j`j], {0, (rmax - r0)/dr + 1}]
avNr = ParallelTable[{{r,N[Mean[Table[Length[Nearest[data, datarmin1[[i]], {Infinity, r}]], {i, 1,Length[datarmin1]}]]]}, j`j++;}[[1]], {r, 0.01, 1.01, 0.05},DistributedContexts -> Automatic] // AbsoluteTiming

Then, one may fit this to a power-law, ie $N(r)\sim r^D$ and determine that indeed a uniform distribution of galaxies has a fractal dimension of D=3:

ff = FindFit[avNr[[2]], a r^D2, {a, D2}, r]

and a plot of the result:

enter image description here

The code can be adjusted for other distributions, as the one given by the OP (which will have a different fractal dimension).

  • $\begingroup$ What is $r$ in your code? Is that the cell size? $\endgroup$
    – Cham
    Commented Aug 5, 2022 at 20:26
  • $\begingroup$ It's the radius of a sphere with the observer at the centre. This kind of calculations are mainly done in spherical coordinates. $\endgroup$
    – Hans Olo
    Commented Aug 5, 2022 at 20:27
  • 1
    $\begingroup$ I'll have to study this, since many of the commands are unknown to me. $\endgroup$
    – Cham
    Commented Aug 5, 2022 at 20:28
  • $\begingroup$ Hans, please, see my edit, in my question. $\endgroup$
    – Cham
    Commented Aug 5, 2022 at 22:42
  • $\begingroup$ I think my answer still stands, as the points are by default deltas. In this notation $r \simeq \ell$. $\endgroup$
    – Hans Olo
    Commented Aug 6, 2022 at 6:37

Here is an approach that I hope gets after your objective of

I would like to define the local number of particles per unit volume V=ℓ3, on a variable scale ℓ, and plot it as level curves (i.e surfaces) in 3D.

(* Generate a universe *)
galaxies[p_, q_, r_] := Module[{pts = RandomReal[{-1, 1}, {p, 3}]}, 
  pts = Nest[Join[#, {RandomChoice[#] + r RandomReal[{-1, 1}, 3]}] &, pts, q];
universe = galaxies[500, 10000, 0.05];

(* Number of cells in each dimension of a n x n x n grid to hold densities *)
n = 101;
(* Set a size of cube for the calculation of density *)
(* Galaxies within s in each dimension are included for the
   calculation of density for each grid cell *)
s = 0.05;
(* Array to hold the densities *)
d = ConstantArray[0, {n, n, n}];
(* Add each galaxy to the cells within s of each galaxy location *)
Do[{x, y, z} = universe[[i]];
 ix0 = Max[1, Ceiling[1/2 (1 + n - (x - s) + n (x - s))]]; 
 ix1 = Min[n, Floor[(1 + n - (x + s) + n (x + s))/2]];
 iy0 = Max[1, Ceiling[1/2 (1 + n - (y - s) + n (y - s))]]; 
 iy1 = Min[n, Floor[(1 + n - (y + s) + n (y + s))/2]];
 iz0 = Max[1, Ceiling[1/2 (1 + n - (z - s) + n (z - s))]]; 
 iz1 = Min[n, Floor[(1 + n - (z + s) + n (z + s))/2]];
    d[[ix, iy, iz]] = d[[ix, iy, iz]] + 1/(2 s)^3, {ix, ix0, ix1}], {iy, iy0, iy1}], {iz, iz0, iz1}],
 {i, Length[universe]}]

(* Create array to hold coordinates and density for cell in d *)
(* This is for the required input format to ListContourPlot3D *)
density = ConstantArray[{0, 0, 0, 0}, n^3];
m = 0;
(* Convert cell indices to locations and add to the density array *)
Do[x = (-1 + 2 i - n)/(-1 + n); 
 Do[y = (-1 + 2 j - n)/(-1 + n); 
  Do[z = (-1 + 2 k - n)/(-1 + n);
   m = m + 1; 
   density[[m]] = {x, y, z, d[[i, j, k]]}, {i, n}], {j, n}], {k, n}]

(* Take a look at the resulting densities to consider *)
Histogram[density[[All, 4]], "FreedmanDiaconis", Frame -> True, 
 FrameLabel -> {"Density", "Count"}]

Histogram of densities

(* Choose a density of interest for contour *)
contourDensity = 5000;

(* ContourPlot3D *)
LBoxRatios -> {1, 1, 1}, MaxPlotPoints -> 100]

3D contour plot of galaxy density

(* As a check on what is produced by ContourPlot3D, 
   show all of the grid points with a density of at least countourDensity *)
h = Select[density, #[[4]] > contourDensity &];
ListPointPlot3D[h[[All, {1, 2, 3}]], BoxRatios -> {1, 1, 1}]

Points on grid where density exceeds 5000


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.