From 5b8aa0c2cb6f187eba957de9400913b83bbf539a Mon Sep 17 00:00:00 2001
From: Arun Isaac
Date: Fri, 19 Mar 2021 21:46:38 +0530
Subject: Factor out random planar angle generation.

* contrib/cone-vector.py (random_planar_angle_cdf,
random_vector_on_disk): New functions.
(random_vector_on_spherical_cap): Rename to
random_vector_on_spherical_cap_cdf. Call random_planar_angle_cdf and
random_vector_on_disk.
(sample_code): Call random_vector_on_spherical_cap_cdf instead of
random_vector_on_spherical_cap.
---
 contrib/cone-vector.py | 22 ++++++++++++++++------
 1 file changed, 16 insertions(+), 6 deletions(-)

diff --git a/contrib/cone-vector.py b/contrib/cone-vector.py
index 8d1877e..65e08ee 100644
--- a/contrib/cone-vector.py
+++ b/contrib/cone-vector.py
@@ -56,24 +56,34 @@ def rotate_from_nth_canonical(x, axis):
             - axisn*(xn*s + a*(axisn - 1))/b
     return x
 
-def random_vector_on_spherical_cap(axis, maximum_planar_angle):
-    """Return a random vector uniformly distributed on a spherical cap."""
+def random_planar_angle_cdf(maximum_planar_angle, dim):
+    """Return a random planar angle using inverse transform sampling."""
+    return solid_angle_fraction2planar_angle(
+        planar_angle2solid_angle_fraction(maximum_planar_angle, dim)*random(),
+        dim)
+
+def random_vector_on_disk(axis, planar_angle):
+    """Return a random vector uniformly distributed on the periphery of a
+disk."""
     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] = sin(planar_angle) * random_vector_on_sphere(dim - 1)
     x[-1] = cos(planar_angle)
     return rotate_from_nth_canonical(x, axis)
 
+def random_vector_on_spherical_cap_cdf(axis, maximum_planar_angle):
+    """Return a random vector uniformly distributed on a spherical
+cap. The random planar angle is generated using inverse transform
+sampling."""
+    return random_vector_on_disk(axis, random_planar_angle_cdf(maximum_planar_angle, axis.size))
+
 def sample_code():
     """Run some sample code testing the defined functions."""
     dim = 100
     maximum_planar_angle = pi/3
     axis = ones(dim)
     axis = axis/norm(axis)
-    print(random_vector_on_spherical_cap(axis, maximum_planar_angle))
+    print(random_vector_on_spherical_cap_cdf(axis, maximum_planar_angle))
 
 if __name__ == '__main__':
     sample_code()
-- 
cgit v1.2.3