I am looking to be able to generate a random uniform sample of particle locations that fall within a spherical volume.
The image below (courtesy of http://nojhan.free.fr/metah/) shows what I am looking for. This is a slice through the sphere, showing a uniform distribution of points:
This is what I am currently getting:
You can see that there is a cluster of points at the center due to the conversion between spherical and Cartesian coordinates.
The code I am using is:
def new_positions_spherical_coordinates(self):
radius = numpy.random.uniform(0.0,1.0, (self.number_of_particles,1))
theta = numpy.random.uniform(0.,1.,(self.number_of_particles,1))*pi
phi = numpy.arccos(1-2*numpy.random.uniform(0.0,1.,(self.number_of_particles,1)))
x = radius * numpy.sin( theta ) * numpy.cos( phi )
y = radius * numpy.sin( theta ) * numpy.sin( phi )
z = radius * numpy.cos( theta )
return (x,y,z)
Below is some MATLAB code that supposedly creates a uniform spherical sample, which is similar to the equation given by http://nojhan.free.fr/metah. I just can't seem to decipher it or understand what they did.
function X = randsphere(m,n,r)
% This function returns an m by n array, X, in which
% each of the m rows has the n Cartesian coordinates
% of a random point uniformly-distributed over the
% interior of an n-dimensional hypersphere with
% radius r and center at the origin. The function
% 'randn' is initially used to generate m sets of n
% random variables with independent multivariate
% normal distribution, with mean 0 and variance 1.
% Then the incomplete gamma function, 'gammainc',
% is used to map these points radially to fit in the
% hypersphere of finite radius r with a uniform % spatial distribution.
% Roger Stafford - 12/23/05
X = randn(m,n);
s2 = sum(X.^2,2);
X = X.*repmat(r*(gammainc(s2/2,n/2).^(1/n))./sqrt(s2),1,n);
I would greatly appreciate any suggestions on generating a truly uniform sample from a spherical volume in Python.
There seem to be plenty of examples showing how to sample from a uniform spherical shell, but that seems to be easier an easier problem. The issue has to do with the scaling - there should be fewer particles at a radius of 0.1 than at a radius of 1.0 to generate a uniform sample from the volume of the sphere.
Edit: Fixed and removed the fact I asked for normally and I meant uniform.
While I prefer the discarding method for spheres, for completeness I offer the exact solution.
In spherical coordinates, taking advantage of the sampling rule:
phi = random(0,2pi)
costheta = random(-1,1)
u = random(0,1)
theta = arccos( costheta )
r = R * cuberoot( u )
now you have a (r, theta, phi)
group which can be transformed to (x, y, z)
in the usual way
x = r * sin( theta) * cos( phi )
y = r * sin( theta) * sin( phi )
z = r * cos( theta )