Deal in solid angle fractions, not absolute solid angles.

* src/extent-sampling.sc (volume-cone): Use
solid-angle-fraction->planar-angle instead of
solid-angle->planar-angle.
* src/nd-random.sc (planar-angle->solid-angle): Rename to
planar-angle->solid-angle-fraction and return a solid angle fraction.
(solid-angle->planar-angle): Rename to
solid-angle-fraction->planar-angle and accept a solid angle fraction.
(hollow-cone-random-vector): Use planar-angle->solid-angle-fraction
instead of planar-angle->solid-angle.
(beta-inc-unnormalized, incomplete-wallis-integral): Delete functions.
* include/nd-random.h (planar_angle_to_solid_angle): Rename to
planar_angle_to_solid_angle_fraction.
(solid_angle_to_planar_angle): Rename to
solid_angle_fraction_to_planar_angle.
* scm/nsmc/wrap.scm (planar-angle->solid-angle): Rename to
planar-angle->solid-angle-fraction.

2 files changed, 27 insertions, 33 deletions

diff --git a/src/extent-sampling.sc b/src/extent-sampling.sc
index 48e6544..73f0e3d 100644
--- a/src/extent-sampling.sc
+++ b/src/extent-sampling.sc
@@ -142,8 +142,8 @@
   (declare volume double)
   (let* ((dimension (unsigned int) (: mean size))
          (vn double (ln-volume-of-ball dimension))
-         (theta-min double (solid-angle->planar-angle (* omega-min (surface-area-of-ball dimension)) dimension))
-         (theta-max double (solid-angle->planar-angle (* omega-max (surface-area-of-ball dimension)) dimension)))
+         (theta-min double (solid-angle-fraction->planar-angle omega-min dimension))
+         (theta-max double (solid-angle-fraction->planar-angle omega-max dimension)))
     (with-rstats stats
       (with-vector x dimension
         (for-i i number-of-samples
diff --git a/src/nd-random.sc b/src/nd-random.sc
index 7044758..79f6cde 100644
--- a/src/nd-random.sc
+++ b/src/nd-random.sc
@@ -60,51 +60,45 @@ sphere. Write the result to X."
                                            (* a (- mun 1))))
                                  b)))))))))
 
-(define (beta-inc-unnormalized a b x) ((static double) double double double)
-  (return (* (gsl-sf-beta-inc a b x)
-             (gsl-sf-beta a b))))
-
-(define (incomplete-wallis-integral theta m) ((static double) double (unsigned int))
-  "Return the incomplete Wallis integral \\int_0^\\theta \\sin^m x
-dx. THETA should be in [0,pi]."
-  (cond
-   ((or (< theta 0) (> theta M-PI))
-    (GSL-ERROR "Incomplete Wallis integral only allows theta in [0,pi]" GSL-EDOM))
-   ((< theta M-PI-2)
-    (return (/ (beta-inc-unnormalized (* 0.5 (+ m 1)) 0.5 (gsl-pow-2 (sin theta)))
-               2)))
-   (else
-    (return (/ (+ (gsl-sf-beta (* 0.5 (+ m 1)) 0.5)
-                  (beta-inc-unnormalized 0.5 (* 0.5 (+ m 1)) (gsl-pow-2 (cos theta))))
-               2)))))
-
-(define (planar-angle->solid-angle planar-angle dimension) (double double (unsigned int))
-  (return (/ (* 2
-                (pow M-PI (* 0.5 (- dimension 1)))
-                (incomplete-wallis-integral planar-angle (- dimension 2)))
-             (gsl-sf-gamma (* 0.5 (- dimension 1))))))
+(define (planar-angle->solid-angle-fraction theta dimension)
+  (double double (unsigned int))
+  (let* ((alpha double (* 0.5 (- dimension 1)))
+         (beta double 0.5)
+         (x double (gsl-pow-2 (sin theta))))
+    (cond
+     ((or (< theta 0) (> theta M-PI))
+      (GSL-ERROR "planar-angle->solid-angle-fraction only allows theta in [0,pi]" GSL-EDOM))
+     ;; gsl-sf-beta-inc underflows when x = 0. Work around that. TODO:
+     ;; Fix it upstream in gsl.
+     ((< (rerror theta M-PI) DBL-EPSILON)
+      (return 1))
+     ((< theta M-PI-2)
+      (return (* 0.5 (gsl-sf-beta-inc alpha beta x))))
+     (else
+      (return (- 1 (* 0.5 (gsl-sf-beta-inc alpha beta x))))))))
 
-(define (solid-angle->planar-angle solid-angle dimension) (double double (unsigned int))
+(define (solid-angle-fraction->planar-angle solid-angle-fraction dimension)
+  (double double (unsigned int))
   (define (f planar-angle params) (double double void*)
-    (return (- (planar-angle->solid-angle planar-angle dimension)
-               solid-angle)))
+    (return (- (planar-angle->solid-angle-fraction planar-angle dimension)
+               solid-angle-fraction)))
 
   (declare gsl-f (struct-variable gsl-function (function (address-of f))))
   ;; TODO: Equality comparisons to floating point values may be
   ;; problematic. Fix it.
   (cond
-   ((= solid-angle 0) (return 0))
-   ((= solid-angle (surface-area-of-ball dimension)) (return M-PI))
+   ((= solid-angle-fraction 0) (return 0))
+   ((= solid-angle-fraction 1) (return M-PI))
    (else (return (bisection (address-of gsl-f) 0 M-PI)))))
 
 (define (hollow-cone-random-vector r mean theta-min theta-max x)
   (void (const gsl-rng*) (const gsl-vector*) double double gsl-vector*)
   ;; Generate random vector around the nth canonical basis vector.
   (let* ((n size-t (: x size))
-         (theta double (solid-angle->planar-angle
+         (theta double (solid-angle-fraction->planar-angle
                         (gsl-ran-flat r
-                                      (planar-angle->solid-angle theta-min n)
-                                      (planar-angle->solid-angle theta-max n))
+                                      (planar-angle->solid-angle-fraction theta-min n)
+                                      (planar-angle->solid-angle-fraction theta-max n))
                         n)))
     (gsl-ran-dir-nd r (- n 1) (: x data))
     (gsl-vector-scale x (sin theta))