diff options
author | Arun Isaac | 2021-03-26 14:43:07 +0530 |
---|---|---|
committer | Arun Isaac | 2021-03-26 14:48:48 +0530 |
commit | b6e2f7875029d1a5a0f6f7a5761b5f3af1f7b448 (patch) | |
tree | f9a470f6e947a0cfc1bec1b542318bd8b66a857b /sambal | |
parent | 05f6d0e94e5603caa769137508cf005bea3c9adc (diff) | |
download | sambal-b6e2f7875029d1a5a0f6f7a5761b5f3af1f7b448.tar.gz sambal-b6e2f7875029d1a5a0f6f7a5761b5f3af1f7b448.tar.lz sambal-b6e2f7875029d1a5a0f6f7a5761b5f3af1f7b448.zip |
Simplify directory structure.
* setup.cfg: Specify the sambal package explicitly.
* src/sambal: Move to sambal.
Diffstat (limited to 'sambal')
-rw-r--r-- | sambal/sambal.py | 68 |
1 files changed, 68 insertions, 0 deletions
diff --git a/sambal/sambal.py b/sambal/sambal.py new file mode 100644 index 0000000..f4d64f4 --- /dev/null +++ b/sambal/sambal.py @@ -0,0 +1,68 @@ +# sambal --- Sample balls, spheres, spherical caps +# Copyright © 2021 Arun I <arunisaac@systemreboot.net> +# Copyright © 2021 Murugesan Venkatapathi <murugesh@iisc.ac.in> +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see +# <https://www.gnu.org/licenses/>. + +from numpy import cos, dot, empty, log, sin, sqrt, pi +from numpy.random import randn, random +from numpy.linalg import norm + +def random_on_sphere(dim): + """Return a random vector uniformly distributed on the unit sphere.""" + x = randn(dim) + return x / norm(x) + +def rotate_from_nth_canonical(x, axis): + """Rotate vector from around the nth canonical axis to the given axis. + + """ + xn = x[-1] + axisn = axis[-1] + if axisn != 1: + b = norm(axis[:-1]) + a = (dot(x, axis) - xn*axisn) / b + s = sqrt(1 - axisn**2) + x = x + (xn*s + a*(axisn - 1))/b * axis + x[-1] = x[-1] + xn*(axisn - 1) - a*s \ + - axisn*(xn*s + a*(axisn - 1))/b + return x + +def random_on_disk(axis, planar_angle): + """Return a random vector uniformly distributed on the periphery of a +disk. + + """ + dim = axis.size + x = empty(dim) + x[:-1] = sin(planar_angle) * random_on_sphere(dim - 1) + x[-1] = cos(planar_angle) + return rotate_from_nth_canonical(x, axis) + +def random_on_cap(axis, maximum_planar_angle): + """Return a random vector uniformly distributed on a spherical +cap. The random planar angle is generated using rejection sampling. + + """ + # 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) |