Dispersing points on the surface of a sphere

Dispersing or choosing points, uniformly or randomly, on the surface of a sphere is a task that someone like me who does molecular simulations, runs into, every once in a while. In some form or shape, this problem is encountered, for example, while modeling a 3D random walk, the motion of a Brownian particle in space, decorating a nanoparticle with stabilizers etc.

While the task is not hard, there are some elegant methods and potential pitfalls that a person doing this for the first time should be aware of.

As I mentioned earlier, one may choose points randomly, or uniformly on the surface (of a unit sphere, here, but can be generalized trivially).

1. Random:

The key point here is that it is incorrect to choose points by selecting angles phi and theta corresponding to spherical coordinates from uniform distributions. This incorrectly concentrates points near the poles.

The correct method is to select two random numbers u and v from a uniform distribution, and construct the two angles as follows:

u = rand();            // uniform random number between

v = rand();            // zero and one.

phi = 2 * pi * u;      // pi is 3.141...

theta = acos(2*v - 1)  // acos is inverse cosine

This is explained with figures here.

2. Uniform:

This is trickier than I thought when I first attempted to do it. Here is an old link (circa 1998, but still good!) which touches upon some of the subtleties. The trouble starts with what "uniform" means, for an arbitrary number of points to be distributed.

If we decide that uniform means a distribution which "maximizes the minimum separation" between n points on the sphere, we can start going somewhere.

There are a number of algorithms that can try to solve this problem, as described on this excellent page. Here is C++ code that I wrote while implementing the golden spiral rule from that page.

//
// Spray n points uniformly on the surface of a sphere 

// of radius "radius"
//

// The positions are returned as an array of Vectors

// where, struct Vector { double x; double y; double z }

//

// Based on algorithm from 

// http://www.xsi-blog.com/archives/115

//

void SprayPointsSphere(int n, double radius, Vector * p)
{
  double inc = PI * (3.0 - sqrt(5.0));
  double off = 2.0/n;

  for(int k = 0; k < n; k++) {
    double y   = k * off - 1 + (off/2);
    double r   = sqrt(1 - y*y);
    double phi = k * inc;

    p[k].x = cos(phi)*r * radius;

    p[k].y = y * radius;

    p[k].z = sin(phi)*r * radius);

  }
}