aboutsummaryrefslogtreecommitdiff
path: root/contrib
diff options
context:
space:
mode:
authorArun Isaac2021-03-15 14:52:03 +0530
committerArun Isaac2021-03-15 15:27:24 +0530
commitd771e8524094493f8a7c8b25b660b18acf837f92 (patch)
tree706db3a717ed1d3cf5a36d5a2591f5af730216d7 /contrib
parentd5ab1ba4a4373821bc27a35f2a5f7dcb9e955dd4 (diff)
downloadnsmc-d771e8524094493f8a7c8b25b660b18acf837f92.tar.gz
nsmc-d771e8524094493f8a7c8b25b660b18acf837f92.tar.lz
nsmc-d771e8524094493f8a7c8b25b660b18acf837f92.zip
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.
Diffstat (limited to 'contrib')
-rw-r--r--contrib/cone-vector.py12
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