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, 57 insertions, 0 deletions

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.