This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: gsl_rand_dir_3d
On Mon, Aug 04, 2003 at 03:27:00PM -0400, Stuart O Anderson wrote:
> After looking at the algorith a while longer I've begun to suspect it has
> a more serious flaw than the s==0 condition. If the Z component is
> uniformly distributed (since P(x*x+y*y < s | s<=1 ) = s) and the final
> vector's orientation in the X-Y plane is also uniform and independent of Z
> then the probability of selecting a final vector with -.5<Z<.5 is equal
> to the probability of selecting a final vector with -1<=Z<=-.5 or .5<=Z<=1.
> Since the surface area of these two sections is not equal but the
> probability of selecting a vector in each section is equal, the
> distribution cannot be uniform over the surface of the sphere, since the
> probability of selecting a vector in any two regions of equal surface
> area must be equal in a uniform distribution.
I haven't studied the method used in the gsl_rand_dir_3d code (so I
should probably keep quiet), but here is the (almost unedited) code
I used in the past to generate a random set of Euler angles. ran()
gives a random number from 0..1 and I use acos() to get the proper
distribution of thetas. For the current application, psi isn't
necessary. It has the advantage that any set of random numbers will
work, while it uses an acos() instead of a sqrt() for the transform.
I remember plotting a set of points they looked pretty evenly
distributed on a sphere. I know I used NR to figure out how to do
it, but even they acknowledge that "copyright does not protect
ideas".
/* orient - Choose a random orientation for the buckyball.
* This is an euler-angle rotation.
* 0 < phi < 2PI , 0 < theta < PI , 0 < psi < 2PI
* phi and psi are chosen with a uniform distribution.
* theta is weighted by sin(theta)
* sin(theta) = dx/dtheta (NR, section 7.2) ->
* x(theta) = -cos(theta) ->
* theta(x) = acos(-x)
*/
int orient(double *phi, double *theta, double *psi)
{
*phi = 2 * PI * ran();
*theta = acos(ran()*2 - 1);
*psi = 2 * PI * ran();
return 0;
}
In evaluating whether these things work, remember that the element
of surface area is
dA = sin(theta) dtheta dphi
This assumes that phi is the "longitude", theta is the "lattitude,
theta=0 is the "north pole", theta=pi/2 is the "equator",
and theta=pi is the "south pole".
--
Jeff Spirko spirko@lehigh.edu spirko@yahoo.com WD3V |=>
The study of non-linear physics is like the study of non-elephant biology.
All theoretical chemistry is really physics;
and all theoretical chemists know it. -- Richard P. Feynman