diff options
| -rw-r--r-- | src/sambal/sambal.py | 22 |
1 files changed, 10 insertions, 12 deletions
diff --git a/src/sambal/sambal.py b/src/sambal/sambal.py index 3a35438..c42ea44 100644 --- a/src/sambal/sambal.py +++ b/src/sambal/sambal.py @@ -39,17 +39,6 @@ def rotate_from_nth_canonical(x, axis): - axisn*(xn*s + a*(axisn - 1))/b return x -def random_planar_angle_pdf(maximum_planar_angle, dim): - """Return a random planar angle using rejection sampling.""" - # We apply the log function just to prevent the floats from - # underflowing. - box_height = (dim-2)*log(sin(min(maximum_planar_angle, pi/2))) - while True: - theta = maximum_planar_angle*random() - f = box_height + log(random()) - if f < (dim-2)*log(sin(theta)): - return theta - def random_on_disk(axis, planar_angle): """Return a random vector uniformly distributed on the periphery of a disk.""" @@ -64,4 +53,13 @@ def random_on_cap(axis, maximum_planar_angle): cap. The random planar angle is generated using rejection sampling. """ - return random_on_disk(axis, random_planar_angle_pdf(maximum_planar_angle, axis.size)) + # We apply the log function just to prevent the floats from + # underflowing. + dim = axis.size + box_height = (dim-2)*log(sin(min(maximum_planar_angle, pi/2))) + while True: + theta = maximum_planar_angle*random() + f = box_height + log(random()) + if f < (dim-2)*log(sin(theta)): + break + return random_on_disk(axis, theta) |