update test-find-peaks to plot cerenkov rings

This commit updates the test-find-peaks script to plot Cerenkov rings for each
of the peaks. It also updates the script to use quad to find the position
instead of using the MC information. Finally, I added a -n argument to the
script to specify how many peaks to draw.

1 files changed, 48 insertions, 0 deletions

diff --git a/src/misc.c b/src/misc.c
index a36e9c9..56be84f 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -256,6 +256,54 @@ double trapz(const double *y, double dx, size_t n)
     return sum*dx/2.0;
 }
 
+/* Compute the first intersection of a ray starting at `pos` with direction
+ * `dir` and a sphere centered at the origin with radius `R`. The distance to
+ * the intersection is stored in `l`.
+ *
+ * Returns 1 if the ray intersects the sphere, and 0 if it doesn't.
+ *
+ * Example:
+ *
+ *     double l;
+ *     double pos[0] = {0,0,0};
+ *     double dir[3] = {1,0,0};
+ *
+ *     if (intersect_sphere(pos,dir,PSUP_RADIUS,&l)) {
+ *         hit[0] = pos[0] + l*dir[0];
+ *         hit[1] = pos[1] + l*dir[1];
+ *         hit[2] = pos[2] + l*dir[2];
+ *         printf("ray intersects sphere at %.2f %.2f %.2f\n", hit[0], hit[1], hit[2]);
+ *     } else {
+ *         printf("ray didn't intersect sphere\n");
+ *     }
+ *
+ */
+int intersect_sphere(double *pos, double *dir, double R, double *l)
+{
+    double b, c;
+
+    b = 2*DOT(dir,pos);
+    c = DOT(pos,pos) - R*R;
+
+    if (b*b - 4*c <= 0) {
+        /* Ray doesn't intersect the sphere. */
+        *l = 0.0;
+        return 0;
+    }
+
+    /* First, check the shorter solution. */
+    *l = (-b - sqrt(b*b - 4*c))/2;
+
+    /* If the shorter solution is less than 0, check the second solution. */
+    if (*l < 0)
+        *l = (-b + sqrt(b*b - 4*c))/2;
+
+    /* If the distance is still negative, we didn't intersect the sphere. */
+    if (*l < 0) return 0;
+
+    return 1;
+}
+
 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