update solid angle test and fix a bug in the solid angle calculation

2 files changed, 9 insertions, 4 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);
     }
 }
diff --git a/test.c b/test.c
index 059fd57..1ede905 100644
--- a/test.c
+++ b/test.c
@@ -52,7 +52,7 @@ int test_solid_angle(char *err)
 
         solid_angle = get_solid_angle(pos,pmt,n,r);
 
-        if (!isclose(solid_angle, solid_angle_results[i].omega, 1e-4, 0)) {
+        if (!isclose(solid_angle, solid_angle_results[i].omega, 1e-5, 0)) {
             sprintf(err, "solid angle = %.2f, but expected %.2f", solid_angle, solid_angle_results[i].omega);
             return 1;
         }