The method described in the article does not produce a spherically symmetric distribution since each value of z is equally likely. This means that a small circle near the north pole will contain as many points as the equator. Instead the different values for z should be selected with a probability that scales like 1/r where r is the radius of that circle (i.e. r = sin(theta)).
Alternatively you can generate random x, y, and z coordinates between -1 and 1 and toss it if it falls outside of the sphere (i.e. if x^2+y^2+z^2 > 1). Then just renormalized the coordinates so that they fall on the sphere.
The benefit of this is that the method easily generalizes to arbitrary dimension.
Alternatively you can generate random x, y, and z coordinates between -1 and 1 and toss it if it falls outside of the sphere (i.e. if x^2+y^2+z^2 > 1). Then just renormalized the coordinates so that they fall on the sphere.
The benefit of this is that the method easily generalizes to arbitrary dimension.