diff options
Diffstat (limited to 'solid_angle.c')
| -rw-r--r-- | solid_angle.c | 11 |
1 files changed, 8 insertions, 3 deletions
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); } } |