update path integral to use a fixed number of points

I noticed when fitting electrons that the cquad integration routine was not
very stable, i.e. it would return different results for *very* small changes in
the fit parameters which would cause the fit to stall.

Since it's very important for the minimizer that the likelihood function not
jump around, I am switching to integrating over the path by just using a fixed
number of points and using the trapezoidal rule. This seems to be a lot more
stable, and as a bonus I was able to combine the three integrals (direct
charge, indirect charge, and time) so that we only have to do a single loop.
This should hopefully make the speed comparable since the cquad routine was
fairly effective at only using as many function evaluations as needed.

Another benefit to this approach is that if needed, it will be easier to port
to a GPU.

1 files changed, 20 insertions, 0 deletions

diff --git a/src/misc.c b/src/misc.c
index ba0a090..d20bc26 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -218,6 +218,26 @@ static struct {
     {100,363.73937555556347},
 };
 
+double trapz(const double *y, double dx, size_t n)
+{
+    /* Returns the integral of `y` using the trapezoidal rule assuming a
+     * constant grid spacing `dx`.
+     *
+     * See https://en.wikipedia.org/wiki/Trapezoidal_rule. */
+    size_t i;
+    double sum = 0.0;
+
+    if (n < 2) return 0.0;
+
+    sum = y[0];
+    for (i = 1; i < n-1; i++) {
+        sum += 2*y[i];
+    }
+    sum += y[n-1];
+
+    return sum*dx/2.0;
+}
+
 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