Implement simplified cone sampling algorithm.

* contrib/cone-vector.py: Don't import tan.
(random_vector_on_spherical_cap): Implement simplified algorithm that
directly samples the surface of the sphere instead of sampling a disk
and projecting it onto the surface.

1 files changed, 6 insertions, 6 deletions

diff --git a/contrib/cone-vector.py b/contrib/cone-vector.py
index 9de71d8..1f28d07 100644
--- a/contrib/cone-vector.py
+++ b/contrib/cone-vector.py
@@ -16,7 +16,7 @@
 # along with this program.  If not, see
 # <https://www.gnu.org/licenses/>.
 
-from numpy import arcsin, cos, dot, empty, ones, sin, sqrt, tan, pi, where, zeros
+from numpy import arcsin, cos, dot, empty, ones, sin, sqrt, pi, where, zeros
 from numpy.random import randn, random
 from numpy.linalg import norm
 from scipy.special import betainc, betaincinv, gamma
@@ -54,12 +54,12 @@ def rotate_from_nth_canonical (x, axis):
 def random_vector_on_spherical_cap (axis, maximum_planar_angle):
     dim = axis.size
     maximum_solid_angle_fraction = planar_angle2solid_angle_fraction(maximum_planar_angle, dim)
+    solid_angle_fraction = maximum_solid_angle_fraction*random()
+    planar_angle = solid_angle_fraction2planar_angle(solid_angle_fraction, dim)
     x = empty(dim)
-    x[:-1] = random_vector_on_sphere(dim - 1) \
-        * cos(maximum_planar_angle) \
-        * tan(solid_angle_fraction2planar_angle(maximum_solid_angle_fraction*random(), dim))
-    x[-1] = cos(maximum_planar_angle)
-    return rotate_from_nth_canonical(x / norm(x), axis)
+    x[:-1] = sin(planar_angle) * random_vector_on_sphere(dim - 1)
+    x[-1] = cos(planar_angle)
+    return rotate_from_nth_canonical(x, axis)
 
 # Sample code exercising the above functions
 dim = 100