aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorArun Isaac2021-03-16 15:01:45 +0530
committerArun Isaac2021-03-16 15:16:37 +0530
commit30104d8ee20733cf4e2195ef9a6cd202a7a257e5 (patch)
treed79c075e5091992a2c57d54d64f84d19b339e3c0 /src
parentb90b54cd3f8d9242ff6d1189e65eaff7212fe44e (diff)
downloadnsmc-30104d8ee20733cf4e2195ef9a6cd202a7a257e5.tar.gz
nsmc-30104d8ee20733cf4e2195ef9a6cd202a7a257e5.tar.lz
nsmc-30104d8ee20733cf4e2195ef9a6cd202a7a257e5.zip
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.
Diffstat (limited to 'src')
-rw-r--r--src/extent-sampling.sc4
-rw-r--r--src/nd-random.sc56
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))