From ac9bc989fea328db6db1820d2676cbbf9335584f Mon Sep 17 00:00:00 2001 From: tlatorre Date: Wed, 4 Jul 2018 14:31:33 -0400 Subject: update solid angle test and fix a bug in the solid angle calculation --- solid_angle.c | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) (limited to 'solid_angle.c') diff --git a/solid_angle.c b/solid_angle.c index f89d65e..d84ff22 100644 --- a/solid_angle.c +++ b/solid_angle.c @@ -23,13 +23,18 @@ double get_solid_angle(double *pos, double *pmt, double *n, double r) Rmax = sqrt(L*L + (r0+r)*(r0+r)); k = sqrt(4*r0*r/(L*L + pow(r0+r,2))); - if (r0 <= r) { + if (fabs(r0-r) < 1e-9) { + /* If r0 is very close to r, the GSL routines below will occasionally + * throw a domain error. */ K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE); - P = gsl_sf_ellint_Pcomp(k, a2, GSL_PREC_DOUBLE); + return M_PI - 2*L*K/Rmax; + } else if (r0 <= r) { + K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE); + P = gsl_sf_ellint_Pcomp(k, -a2, GSL_PREC_DOUBLE); return 2*M_PI - 2*L*K/Rmax + (2*L/Rmax)*(r0-r)*P/(r0+r); } else { K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE); - P = gsl_sf_ellint_Pcomp(k, a2, GSL_PREC_DOUBLE); + P = gsl_sf_ellint_Pcomp(k, -a2, GSL_PREC_DOUBLE); return -2*L*K/Rmax + (2*L/Rmax)*(r0-r)*P/(r0+r); } } -- cgit