diff options
| -rw-r--r-- | src/Makefile | 2 | ||||
| -rw-r--r-- | src/misc.c | 57 | ||||
| -rw-r--r-- | src/misc.h | 1 | ||||
| -rw-r--r-- | src/test.c | 137 |
4 files changed, 196 insertions, 1 deletions
diff --git a/src/Makefile b/src/Makefile index 3f04dd8..e15b07d 100644 --- a/src/Makefile +++ b/src/Makefile @@ -20,7 +20,7 @@ test-likelihood: test-likelihood.o muon.o random.o optics.o quantum_efficiency.o test-path: test-path.o mt19937ar.o vector.o path.o random.o misc.o -test-charge: test-charge.o sno_charge.o misc.o +test-charge: test-charge.o sno_charge.o misc.o vector.o test-zebra: test-zebra.o zebra.o pack2b.o @@ -2,6 +2,7 @@ #include <math.h> #include <stdlib.h> /* for size_t */ #include <gsl/gsl_sf_gamma.h> +#include "vector.h" static struct { int n; @@ -217,6 +218,62 @@ static struct { {100,363.73937555556347}, }; +void get_path_length(double *pos1, double *pos2, double R, double *l1, double *l2) +{ + /* Returns the path length inside and outside a circle of radius `R` for a + * ray starting at position `pos1` and ending at position `pos2`. + * + * The path length inside the sphere is stored in `l1` and the path length + * outside the sphere is stored in `l2`. */ + double dir[3], l, b, c, d1, d2; + + /* Calculate the vector from `pos1` to `pos2`. */ + SUB(dir,pos2,pos1); + + l = NORM(dir); + + normalize(dir); + + b = 2*DOT(dir,pos1); + c = DOT(pos1,pos1) - R*R; + + if (b*b - 4*c <= 0) { + /* Ray doesn't intersect the sphere. */ + *l1 = 0.0; + *l2 = l; + return; + } + + d1 = (-b + sqrt(b*b - 4*c))/2; + d2 = (-b - sqrt(b*b - 4*c))/2; + + if (d1 < 0) { + /* Ray also doesn't intersect sphere. */ + *l1 = 0.0; + *l2 = l; + } else if (d1 >= l && d2 < 0) { + /* Ray also doesn't intersect sphere. */ + *l1 = l; + *l2 = 0.0; + } else if (d2 < 0) { + /* Ray intersects sphere once. */ + *l1 = d1; + *l2 = l-d1; + } else if (d1 >= l && d2 >= l) { + /* Ray doesn't intersect the sphere. */ + *l1 = 0.0; + *l2 = l; + } else if (d1 >= l && d2 < l) { + /* Ray intersects the sphere once. */ + *l2 = d1; + *l1 = l-d1; + } else if (d1 < l && d2 < l) { + /* Ray intersects the sphere twice. */ + *l1 = d1-d2; + *l2 = l-(d1-d2); + } +} + double ln(unsigned int n) { /* Returns the logarithm of n. @@ -6,6 +6,7 @@ #define LN_MAX 100 #define LNFACT_MAX 100 +void get_path_length(double *pos1, double *pos2, double R, double *l1, double *l2); double ln(unsigned int n); double lnfact(unsigned int n); double kahan_sum(double *x, size_t n); @@ -12,6 +12,7 @@ #include "random.h" #include "pdg.h" #include <gsl/gsl_spline.h> +#include "vector.h" typedef int testFunction(char *err); @@ -511,6 +512,141 @@ err: return 1; } +int test_get_path_length(char *err) +{ + /* Tests the get_path_length() function. */ + size_t i; + double pos1[3], pos2[3], tmp[3], r; + double l1, l2, l1_expected, l2_expected; + + init_genrand(0); + + /* First we test two points within the sphere. */ + + for (i = 1; i < 1000; i++) { + rand_sphere(pos1); + r = genrand_real2(); + MUL(pos1,r); + rand_sphere(pos2); + r = genrand_real2(); + MUL(pos2,r); + + SUB(tmp,pos2,pos1); + + l1_expected = NORM(tmp); + l2_expected = 0.0; + + get_path_length(pos1,pos2,1.0,&l1,&l2); + + if (!isclose(l1, l1_expected, 1e-9, 1e-9)) { + sprintf(err, "get_path_length() returned %.5g for l1, but expected %.5g", l1, l1_expected); + goto err; + } + + if (!isclose(l2, l2_expected, 1e-9, 1e-9)) { + sprintf(err, "get_path_length() returned %.5g for l2, but expected %.5g", l2, l2_expected); + goto err; + } + } + + /* Now we test one point at the origin and the other outside of the sphere. */ + + pos1[0] = 0.0; + pos1[1] = 0.0; + pos1[2] = 0.0; + + for (i = 1; i < 1000; i++) { + rand_sphere(pos2); + r = genrand_real2() + 1.0; + MUL(pos2,r); + + SUB(tmp,pos2,pos1); + + l1_expected = 1.0; + l2_expected = NORM(tmp) - 1.0; + + get_path_length(pos1,pos2,1.0,&l1,&l2); + + if (!isclose(l1, l1_expected, 1e-9, 1e-9)) { + sprintf(err, "get_path_length() returned %.5g for l1, but expected %.5g", l1, l1_expected); + goto err; + } + + if (!isclose(l2, l2_expected, 1e-9, 1e-9)) { + sprintf(err, "get_path_length() returned %.5g for l2, but expected %.5g", l2, l2_expected); + goto err; + } + } + + /* Now we test both points outside the sphere such that the ray doesn't + * intersect the sphere. */ + + for (i = 1; i < 1000; i++) { + pos1[0] = genrand_real2()-0.5; + pos1[1] = genrand_real2()-0.5; + pos1[2] = 1.0 + genrand_real2(); + pos2[0] = genrand_real2()-0.5; + pos2[1] = genrand_real2()-0.5; + pos2[2] = 1.0 + genrand_real2(); + + SUB(tmp,pos2,pos1); + + l1_expected = 0.0; + l2_expected = NORM(tmp); + + get_path_length(pos1,pos2,1.0,&l1,&l2); + + if (!isclose(l1, l1_expected, 1e-9, 1e-9)) { + sprintf(err, "get_path_length() returned %.5g for l1, but expected %.5g", l1, l1_expected); + goto err; + } + + if (!isclose(l2, l2_expected, 1e-9, 1e-9)) { + sprintf(err, "get_path_length() returned %.5g for l2, but expected %.5g", l2, l2_expected); + goto err; + } + } + + /* Now we test both points outside the sphere such that the ray intersects + * the sphere. */ + + for (i = 1; i < 1000; i++) { + /* Pick a random point outside the sphere. */ + rand_sphere(pos1); + r = genrand_real2() + 1.0; + MUL(pos1,r); + + COPY(pos2,pos1); + MUL(pos2,-1.0); + normalize(pos2); + + r = genrand_real2() + 1.0; + MUL(pos2,r); + + SUB(tmp,pos2,pos1); + + l1_expected = 2.0; + l2_expected = NORM(tmp) - 2.0; + + get_path_length(pos1,pos2,1.0,&l1,&l2); + + if (!isclose(l1, l1_expected, 1e-9, 1e-9)) { + sprintf(err, "get_path_length() returned %.5g for l1, but expected %.5g", l1, l1_expected); + goto err; + } + + if (!isclose(l2, l2_expected, 1e-9, 1e-9)) { + sprintf(err, "get_path_length() returned %.5g for l2, but expected %.5g", l2, l2_expected); + goto err; + } + } + + return 0; + +err: + return 1; +} + struct tests { testFunction *test; char *name; @@ -528,6 +664,7 @@ struct tests { {test_kahan_sum, "test_kahan_sum"}, {test_lnfact, "test_lnfact"}, {test_ln, "test_ln"}, + {test_get_path_length, "test_get_path_length"}, }; int main(int argc, char **argv) |