update test

3 files changed, 39 insertions, 1 deletions

diff --git a/TODO b/TODO
index f8c9736..8efa652 100644
--- a/TODO
+++ b/TODO
@@ -5,4 +5,6 @@
 - write talk for Monday
 - write a section on atmospheric backgrounds in paper
 - check index of refraction calculation against tables in the paper
-- write a test for the solid angle calculation
+- write more tests for the solid angle calculation
+    - use table from paper
+- make sure normal direction is normalized
diff --git a/solid_angle.h b/solid_angle.h
index 99fda37..f4cf96b 100644
--- a/solid_angle.h
+++ b/solid_angle.h
@@ -1,3 +1,8 @@
+#ifndef SOLID_ANGLE_H
+#define SOLID_ANGLE_H
+
 #include <math.h>
 
 double get_solid_angle(double *pos, double *pmt, double *n, double r);
+
+#endif
diff --git a/test.c b/test.c
index 62529ea..059fd57 100644
--- a/test.c
+++ b/test.c
@@ -2,6 +2,24 @@
 #include <math.h>
 #include <stdio.h>
 
+struct solid_angle_results {
+    double L;
+    double r0;
+    double omega;
+} solid_angle_results[] = {
+    {0.5,0.0,3.4732594},
+    {0.5,0.2,3.4184435},
+    {0.5,0.4,3.2435434},
+    {0.5,0.6,2.9185178},
+    {0.5,0.8,2.4122535},
+    {0.5,1.0,1.7687239},
+    {0.5,1.2,1.1661307},
+    {0.5,1.4,0.7428889},
+    {0.5,1.6,0.4841273},
+    {0.5,1.8,0.3287007},
+    {0.5,2.0,0.2324189}
+};
+
 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
@@ -14,6 +32,7 @@ int isclose(double a, double b, double rel_tol, double abs_tol)
 int test_solid_angle(char *err)
 {
     /* Tests the get_solid_angle() function. */
+    int i;
     double pmt[3] = {0,0,0};
     double pos[3] = {0,0,1};
     double n[3] = {0,0,1};
@@ -27,6 +46,18 @@ int test_solid_angle(char *err)
         return 1;
     }
 
+    for (i = 0; i < sizeof(solid_angle_results)/sizeof(struct solid_angle_results); i++) {
+        pos[0] = solid_angle_results[i].r0*r;
+        pos[2] = solid_angle_results[i].L*r;
+
+        solid_angle = get_solid_angle(pos,pmt,n,r);
+
+        if (!isclose(solid_angle, solid_angle_results[i].omega, 1e-4, 0)) {
+            sprintf(err, "solid angle = %.2f, but expected %.2f", solid_angle, solid_angle_results[i].omega);
+            return 1;
+        }
+    }
+
     return 0;
 }