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]

gsl_rand_dir_3d


Hi - I was looking for code to generate uniform spherical distributions
and found gsl_rand_dir_3d().  I think there may be an insignificant error
in the code.  The do while loop that generates the point within a unit
disk rejects the point (0,0), which would be correct if we were trying
to generate a random unit vector in two dimensions (since the next step
would be to normalize the vector to the selected point), but in this case
it is not the correct behavior.  If we can't generate a point a (0,0) the
z component will never get the value -1, and the vector (0,0,-1) can never
be generated.  Allowing (0,0) doesn't adversely affect the rest of the
function.

Is my reasoning correct?  I've included a copy of the function below for
easy reference.

Stuart

void
gsl_ran_dir_3d (const gsl_rng * r, double *x, double *y, double *z)
{
  double s, a;

  /* This is a variant of the algorithm for computing a random point
   * on the unit sphere; the algorithm is suggested in Knuth, v2,
   * 3rd ed, p136; and attributed to Robert E Knop, CACM, 13 (1970),
   * 326.
   */

  /* Begin with the polar method for getting x,y inside a unit circle
   */
  do
    {
      *x = -1 + 2 * gsl_rng_uniform (r);
      *y = -1 + 2 * gsl_rng_uniform (r);
      s = (*x) * (*x) + (*y) * (*y);
    }
  while (s > 1.0 || s == 0.0);

  *z = -1 + 2 * s;		/* z uniformly distributed from -1 to 1 */
  a = 2 * sqrt (1 - s);		/* factor to adjust x,y so that x^2+y^2
				 * is equal to 1-z^2 */
  *x *= a;
  *y *= a;
}


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