add get_path_length()

This commit adds a function called get_path_length() which computes the path
length inside and outside a sphere for a line segment between two points. This
will be useful for calculating the photon absorption for paths which cross the
AV and for computing the time of flight of photons from a track to a PMT.

1 files changed, 137 insertions, 0 deletions

diff --git a/src/test.c b/src/test.c
index 5d813c4..5f9c655 100644
--- a/src/test.c
+++ b/src/test.c
@@ -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)