Disribute particles uniformly inside a sphere with huge radius (10^9) - Matlab

Question

I would like some help with the following:

   L2=10^9;
   nP=10^6;

   x(:,1)=L2.*rand(nP,1);
   y(:,1)=L2.*rand(nP,1);   % Select  initial partice position randomly (only positive values) NS
   z(:,1)= L2.*rand(nP,1);


   for ip=1:nP
      r1(ip,1)= sqrt(x(ip,1)^2 + y(ip,1)^2 + z(ip,1)^2);   % distance from (0,0,0) for each iteration

      while r1(ip,1)>=L2

         x(ip,1)=L2.*rand(1,1);
         y(ip,1)=L2.*rand(1,1);     % Select  initial partice position randomly (only positive values) NS
         z(ip,1)=L2.*rand(1,1);
         r1(ip,1)= sqrt(x(ip,1)^2 + y(ip,1)^2 + z(ip,1)^2);
      end
   end

The problems I face are:

1) Since the radius L2 is very big, when I multiply with L2 with rand I get values very close to the boundaries of the box (10^8 for example) and I want the particles uniformly distributed all around the sphere.

2) Also after I run the while command the resulting distributions are no more uniform.

Why do you say the points are close to the boundaries? And why do you say they are not randomly distributed? Your code looks correct to me. I think you didn’t evaluate the result correctly (except for what I assume are typos, as there are some syntax errors). — Cris Luengo, Mar 10 '20 at 02:29
I corrected the two syntax errors in your code, please make sure it matches the code you were running. — Cris Luengo, Mar 10 '20 at 02:35
If you have a sphere, then you might want to generate your random points using spherical coordinates. Given the shape of a sphere, there is a lot more volume near the outer boundary than near the centre. Some bias towards the boundary is to be expected. — rossum, Mar 10 '20 at 08:59
This question was asking about points distributed on the sphere _surface_ but you may find some inspiration: https://stackoverflow.com/questions/30158433/how-to-generate-random-two-concentric-spheres-synthetic-data/30174449#30174449 — Hoki, Mar 10 '20 at 09:42

pjs · Accepted Answer · 2020-03-11T04:01:23.250

Conceptually this can be accomplished in three stages:

Generate points on the surface of a unit sphere;
Scale each point to fill the interior of the sphere uniformly; and
Rescale them to the desired radius.

In practice the steps can be combined.

Note that I'm not a Matlab user, nor do I have a copy of it, so please treat what follows as pseudocode.

The simplest and most general way to generate points randomly distributed on the surface of a k-dimensional sphere, while avoiding the acceptance/rejection approach, is to generate k independent standard normals and normalize each resulting point by its distance from the origin. For three dimensions:

x = randn
y = randn
z = randn
scale_factor = 1 / sqrt(x^2 + y^2 + z^2)
x = x * scale_factor
y = y * scale_factor
z = z * scale_factor

To scale the points to the interior of the sphere, you want to apply a scale factor between 0 and 1. However, if you scale the factor uniformly the result will not be three dimensionally uniform. Note that a sphere with radius 1/2 has only 1/8 the volume of a unit sphere, so we need to adjust the scale factor accordingly to make the expected number of points proportional to the volume. The correct scale factor for a k-dimensional sphere is to take the k-th root of a uniform(0,1) random number (see Wolfram for a 2-d illustration of this principle), so in your case you want to adjust the scale factor using a cube root. Change the scale_factor above to:

scale_factor = (rand^(1/3)) / sqrt(x^2 + y^2 + z^2)

Finally, to rescale to your actual desired radius of L2, what you actually want is:

scale_factor = L2 * (rand^(1/3)) / sqrt(x^2 + y^2 + z^2)

Disribute particles uniformly inside a sphere with huge radius (10^9) - Matlab

1 Answers1