diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rw-r--r-- | Makefile | 2 | ||||
| -rw-r--r-- | solid_angle.c | 2 | ||||
| -rw-r--r-- | test.c | 44 |
4 files changed, 48 insertions, 1 deletions
@@ -1,3 +1,4 @@ *.swp calculate_limits *.o +test @@ -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); @@ -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; + } +} |