Re: gsl_rand_dir_3d

[James Theiler <jt at lanl dot gov>]

Thanks Stuart,
  I agree with your reasoning.
  Just for grins, I went back to the original citation (Knop, 1970),
and in that algorithm, there is no s==0 condition, just as you
suggest.  Also, for some reason, that algorithm suggests
  a = sqrt((1-(*z)*(*z))/s)
instead of the equivalent a=2*sqrt(1-s) that we have.  Usually, these
old papers have a good numerical reason for these kinds of choices,
but in this case it seems like it would just lead to problems when
s=0.
  Bottom line, I suggest we follow Stuart's recommendation, and
remove the s==0 condition from the do while.

jt

On Sun, 3 Aug 2003, Stuart O Anderson wrote:

] 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;
] }
]

---------------------------------------------
James Theiler                     jt@lanl.gov
MS-B244, NIS-2, LANL        tel: 505/665-5682
Los Alamos, NM 87545        fax: 505/665-4414
----- Space and Remote Sensing Sciences -----