Generate a random sample of points distributed on the surface of a unit sphere

Based on the last approach on this page, you can simply generate a vector consisting of independent samples from three standard normal distributions, then normalize the vector such that its magnitude is 1:

import numpy as npdef sample_spherical(npoints, ndim=3):    vec = np.random.randn(ndim, npoints)    vec /= np.linalg.norm(vec, axis=0)    return vec

For example:

from matplotlib import pyplot as pltfrom mpl_toolkits.mplot3d import axes3dphi = np.linspace(0, np.pi, 20)theta = np.linspace(0, 2 * np.pi, 40)x = np.outer(np.sin(theta), np.cos(phi))y = np.outer(np.sin(theta), np.sin(phi))z = np.outer(np.cos(theta), np.ones_like(phi))xi, yi, zi = sample_spherical(100)fig, ax = plt.subplots(1, 1, subplot_kw={'projection':'3d', 'aspect':'equal'})ax.plot_wireframe(x, y, z, color='k', rstride=1, cstride=1)ax.scatter(xi, yi, zi, s=100, c='r', zorder=10)

The same method also generalizes to picking uniformly distributed points on the unit circle (ndim=2) or on the surfaces of higher-dimensional unit hyperspheres.

Points on the surface of a sphere can be expressed using two spherical coordinates, theta and phi, with 0 < theta < 2pi and 0 < phi < pi.

Conversion formula into cartesian x, y, z coordinates:

x = r * cos(theta) * sin(phi)y = r * sin(theta) * sin(phi)z = r * cos(phi)

where r is the radius of the sphere.

So the program could randomly sample theta and phi in their ranges, at uniform distribution, and generate the cartesian coordinates from it.

But then the points get distributed more densley on the poles of the sphere. In order for points to get uniformly distributed on the sphere surface, phi needs to be chosen as phi = acos(a) where -1 < a < 1 is chosen on an uniform distribution.

For the Numpy code it would be the same as in Sampling uniformly distributed random points inside a spherical volume , except that the variable radius has a fixed value.

Another way that depending on the hardware could be much faster.

Choose a, b, c to be three random numbers each between -1 and 1

Calculate r2 = a^2 + b^2 + c^2

If r2 > 1.0 (=the point isn't in the sphere) or r2 < 0.00001 (=the point is too close to the center, we'll have division by zero while projecting to the surface of the sphere) you discard the values, and pick another set of random a, b, c

Otherwise, you’ve got your random point (relative to center of the sphere):

ir = R / sqrt(r2)x = a * iry = b * irz = c * ir