diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/extent-sampling.sc | 4 | ||||
| -rw-r--r-- | src/nd-random.sc | 56 |
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)) |