diff options
| -rw-r--r-- | contrib/cone-vector.py | 22 |
1 files 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() |