This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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 


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]