From f60a54c593411803a8617c9dcda2f667095d2487 Mon Sep 17 00:00:00 2001
From: Arun Isaac
Date: Wed, 30 Jun 2021 14:42:24 +0530
Subject: Support sampling in sphere.

* src/nd-random.sc (random-vector-in-sphere): New function.
* include/nd-random.h (random_vector_in_sphere): Declare it.
---
 include/nd-random.h | 1 +
 src/nd-random.sc    | 9 +++++++++
 2 files changed, 10 insertions(+)

diff --git a/include/nd-random.h b/include/nd-random.h
index 3d3c298..a260bf1 100644
--- a/include/nd-random.h
+++ b/include/nd-random.h
@@ -27,6 +27,7 @@
 #include <gsl/gsl_vector.h>
 
 void random_direction_vector (const gsl_rng* r, gsl_vector* x);
+void random_vector_in_sphere (const gsl_rng* r, double radius, gsl_vector* x);
 void cone_random_vector (const gsl_rng* r, const gsl_vector* mean, double theta_max, gsl_vector* x);
 void hollow_cone_random_vector (const gsl_rng* r, const gsl_vector* mean, double theta_min, double theta_max, gsl_vector* x);
 
diff --git a/src/nd-random.sc b/src/nd-random.sc
index 79f6cde..69bb93a 100644
--- a/src/nd-random.sc
+++ b/src/nd-random.sc
@@ -30,6 +30,15 @@
 sphere. Write the result to X."
   (gsl-ran-dir-nd r (: x size) (: x data)))
 
+(define (random-vector-in-sphere rng radius x)
+  (void (const gsl-rng*) double gsl-vector*)
+  "Generate a random vector distributed uniformly in a sphere of
+RADIUS. Store the generated vector in X."
+  (random-direction-vector rng x)
+  (gsl-vector-scale x (* radius
+                         (pow (gsl-ran-flat rng 0 1)
+                              (/ (: x size))))))
+
 (define (rotate-from-nth-canonical x orient)
   ((static void) gsl-vector* (const gsl-vector*))
   (let* ((n size-t (: x size))
-- 
cgit v1.2.3