From 30104d8ee20733cf4e2195ef9a6cd202a7a257e5 Mon Sep 17 00:00:00 2001 From: Arun Isaac Date: Tue, 16 Mar 2021 15:01:45 +0530 Subject: 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. --- include/nd-random.h | 4 ++-- scm/nsmc/wrap.scm | 4 ++-- src/extent-sampling.sc | 4 ++-- src/nd-random.sc | 56 ++++++++++++++++++++++---------------------------- 4 files changed, 31 insertions(+), 37 deletions(-) diff --git a/include/nd-random.h b/include/nd-random.h index 9755e57..3d3c298 100644 --- a/include/nd-random.h +++ b/include/nd-random.h @@ -30,7 +30,7 @@ void random_direction_vector (const gsl_rng* r, 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); -double planar_angle_to_solid_angle (double planar_angle, unsigned int dimension); -double solid_angle_to_planar_angle (double solid_angle, unsigned int dimension); +double planar_angle_to_solid_angle_fraction (double planar_angle, unsigned int dimension); +double solid_angle_fraction_to_planar_angle (double solid_angle, unsigned int dimension); #endif diff --git a/scm/nsmc/wrap.scm b/scm/nsmc/wrap.scm index 6ae4089..3e94192 100644 --- a/scm/nsmc/wrap.scm +++ b/scm/nsmc/wrap.scm @@ -64,9 +64,9 @@ (dynamic-func "angle_between_vectors" lib-nsmc) (list '* '*))) -(define-public planar-angle->solid-angle +(define-public planar-angle->solid-angle-fraction (pointer->procedure double - (dynamic-func "planar_angle_to_solid_angle" lib-nsmc) + (dynamic-func "planar_angle_to_solid_angle_fraction" lib-nsmc) (list double int))) ;; Random state 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)) -- cgit v1.2.3