add a test of the get_solid_angle function

4 files changed, 48 insertions, 1 deletions

diff --git a/.gitignore b/.gitignore
index 9e63049..7d3fd46 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,4 @@
 *.swp
 calculate_limits
 *.o
+test
diff --git a/Makefile b/Makefile
index 6148848..98f700d 100644
--- a/Makefile
+++ b/Makefile
@@ -5,6 +5,8 @@ calculate_limits: calculate_limits.c
 
 solid_angle.o: solid_angle.c
 
+test: test.c solid_angle.c
+
 clean:
 	rm calculate_limits
 
diff --git a/solid_angle.c b/solid_angle.c
index 5fc6033..cbc71f1 100644
--- a/solid_angle.c
+++ b/solid_angle.c
@@ -15,7 +15,7 @@ double get_solid_angle(double *pos, double *pmt, double *n, double r)
     dir[1] = pos[1] - pmt[1];
     dir[2] = pos[2] - pmt[2];
 
-    L = dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2];
+    L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]);
     R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
     r0 = sqrt(R*R - L*L);
 
diff --git a/test.c b/test.c
new file mode 100644
index 0000000..62529ea
--- /dev/null
+++ b/test.c
@@ -0,0 +1,44 @@
+#include "solid_angle.h"
+#include <math.h>
+#include <stdio.h>
+
+int isclose(double a, double b, double rel_tol, double abs_tol)
+{
+    /* Returns 1 if a and b are "close". This algorithm is taken from Python's
+     * math.isclose() function.
+     *
+     * See https://www.python.org/dev/peps/pep-0485/. */
+    return fabs(a-b) <= fmax(rel_tol*fmax(fabs(a),fabs(b)),abs_tol);
+}
+
+int test_solid_angle(char *err)
+{
+    /* Tests the get_solid_angle() function. */
+    double pmt[3] = {0,0,0};
+    double pos[3] = {0,0,1};
+    double n[3] = {0,0,1};
+    double r = 1.0;
+    double solid_angle;
+
+    solid_angle = get_solid_angle(pos,pmt,n,r);
+
+    if (!isclose(solid_angle, 2*M_PI*(1-1/sqrt(2)), 1e-9, 0)) {
+        sprintf(err, "solid angle = %.2f, but expected %.2f", solid_angle, 2*M_PI*(1-1/sqrt(2)));
+        return 1;
+    }
+
+    return 0;
+}
+
+int main(int argc, char **argv)
+{
+    char err[256];
+
+    if (!test_solid_angle(err)) {
+        printf("[\033[92mok\033[0m] test_solid_angle\n");
+        return 0;
+    } else {
+        printf("[\033[91mfail\033[0m] test_solid_angle: %s\n", err);
+        return 1;
+    }
+}