a bunch of small changes to speed things up

1 files changed, 28 insertions, 0 deletions

diff --git a/src/sno_charge.c b/src/sno_charge.c
index 1f39706..a019134 100644
--- a/src/sno_charge.c
+++ b/src/sno_charge.c
@@ -14,6 +14,8 @@ double qhi = 100.0;
 size_t nq = 10000;
 static gsl_spline *splines[MAX_PE];
 static gsl_interp_accel *acc;
+static gsl_spline *log_splines[MAX_PE];
+static gsl_interp_accel *log_acc;
 
 static double qmean, qstd;
 
@@ -84,6 +86,23 @@ double spe_pol2exp(double q)
     return (gmratio*(NG1*funpol1 + NG2*funpol2) + NEXP*exp(-q/Q0))/norm;
 }
 
+double log_pq(double q, int n)
+{
+    if (!initialized) {
+        fprintf(stderr, "charge interpolation hasn't been initialized!\n");
+        exit(1);
+    }
+
+    if (n > MAX_PE) {
+        /* Assume the distribution is gaussian by the central limit theorem. */
+        return log(norm(q,n*qmean,sqrt(n)*qstd));
+    }
+
+    if (q < qlo || q > qhi) return -INFINITY;
+
+    return interp1d(q,log_splines[n-1]->x,log_splines[n-1]->y,log_splines[n-1]->size);
+}
+
 double pq(double q, int n)
 {
     if (!initialized) {
@@ -240,6 +259,15 @@ void init_charge(void)
     }
 
     gsl_integration_cquad_workspace_free(w);
+
+    for (i = 1; i <= MAX_PE; i++) {
+        for (j = 0; j < nq; j++) {
+            y[j] = log(splines[i-1]->y[j]);
+        }
+        log_splines[i-1] = gsl_spline_alloc(gsl_interp_linear,nq);
+        gsl_spline_init(log_splines[i-1],x,y,nq);
+    }
+
     free(x);
     free(y);
 }