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.

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
 
diff --git a/src/misc.c b/src/misc.c
index f0e37e3..00f2f6c 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -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.
diff --git a/src/misc.h b/src/misc.h
index 6e9095a..f44125a 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -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);
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)