speed up likelihood function by not calling trapz()

This commit speeds up the likelihood function by integrating the charge along
the track inline instead of creating an array and then calling trapz(). It also
introduces two global variables avg_index_d2o and avg_index_h2o which are the
average indices of refraction for D2O and H2O weighted by the PMT quantum
efficiency and the Cerenkov spectrum.

3 files changed, 140 insertions, 91 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index d29be5e..0228bd8 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -208,31 +208,34 @@ void particle_free(particle *p)
     free(p);
 }
 
-static double get_expected_charge_shower(particle *p, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o, double *q_delta_ray, double *q_indirect_delta_ray)
+static void get_expected_charge_shower(particle *p, double *pos, double *dir, int pmt, double *q, double *reflected, double *q_delta_ray, double *q_indirect_delta_ray, double *t)
 {
-    double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, charge, constant, prob_abs, prob_sct;
+    double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, charge, constant, prob_abs, prob_sct, l_d2o, l_h2o;
 
-    SUB(pmt_dir,pmt_pos,pos);
+    SUB(pmt_dir,pmts[pmt].pos,pos);
 
     normalize(pmt_dir);
 
-    cos_theta_pmt = DOT(pmt_dir,pmt_normal);
+    cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal);
 
-    charge = 0.0;
+    *t = 0.0;
+    *q = 0.0;
+    *reflected = 0.0;
     *q_delta_ray = 0.0;
     *q_indirect_delta_ray = 0.0;
-    *reflected = 0.0;
-    if (cos_theta_pmt > 0) return 0.0;
+    if (cos_theta_pmt > 0) return;
 
     /* Calculate the cosine of the angle between the track direction and the
      * vector to the PMT. */
     cos_theta = DOT(dir,pmt_dir);
 
-    omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r);
+    omega = get_solid_angle_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS);
 
     f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
     f = get_weighted_pmt_response(-cos_theta_pmt);
 
+    get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);
+
     prob_abs = 1.0 - get_fabs_d2o(l_d2o)*get_fabs_h2o(l_h2o)*get_fabs_acrylic(AV_THICKNESS);
     prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o);
 
@@ -244,34 +247,39 @@ static double get_expected_charge_shower(particle *p, double *pos, double *dir,
      * showers from lower energy electrons or for delta rays from lower energy
      * muons. It seems good enough, but in the future it would be nice to
      * parameterize this. */
+    charge = 0.0;
     if (p->shower_photons > 0)
-        charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/n_d2o);
+        charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/avg_index_d2o);
 
     if (p->delta_ray_photons > 0)
-        *q_delta_ray = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/n_d2o);
+        *q_delta_ray = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/avg_index_d2o);
 
     *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
     *q_indirect_delta_ray = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*(*q_delta_ray) + prob_sct*(*q_delta_ray);
 
     *q_delta_ray *= (1.0-prob_abs)*(1.0-prob_sct)*f;
 
-    return (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
+    *t = (l_d2o*avg_index_d2o + l_h2o*avg_index_h2o)/SPEED_OF_LIGHT;
+
+    *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
 }
 
-static double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o)
+static void get_expected_charge(double beta, double theta0, double *pos, double *dir, int pmt, double *q, double *reflected, double *t)
 {
-    double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, charge, prob_abs, prob_sct;
+    double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o;
 
     z = 1.0;
 
     R = NORM(pos);
 
-    n = (R <= AV_RADIUS) ?  n_d2o : n_h2o;
+    n = (R <= AV_RADIUS) ? avg_index_d2o : avg_index_h2o;
 
+    *q = 0.0;
+    *t = 0.0;
     *reflected = 0.0;
-    if (beta < 1/n) return 0.0;
+    if (beta < 1/n) return;
 
-    SUB(pmt_dir,pmt_pos,pos);
+    SUB(pmt_dir,pmts[pmt].pos,pos);
 
     normalize(pmt_dir);
 
@@ -279,18 +287,21 @@ static double get_expected_charge(double x, double beta, double theta0, double *
      * vector to the PMT. */
     cos_theta = DOT(dir,pmt_dir);
 
-    *reflected = 0.0;
-    if (fabs(cos_theta-1.0/(n_d2o*beta))/theta0 > 5) return 0.0;
+    if (fabs(cos_theta-1.0/(n*beta))/theta0 > 5) return;
 
-    cos_theta_pmt = DOT(pmt_dir,pmt_normal);
+    cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal);
 
-    if (cos_theta_pmt > 0) return 0.0;
+    if (cos_theta_pmt > 0) return;
 
-    omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r);
+    omega = get_solid_angle_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS);
 
     f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
     f = get_weighted_pmt_response(-cos_theta_pmt);
 
+    get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);
+
+    *t = (l_d2o*avg_index_d2o + l_h2o*avg_index_h2o)/SPEED_OF_LIGHT;
+
     /* Probability that a photon is absorbed. We calculate this by computing:
      *
      *     1.0 - P(not absorbed in D2O)*P(not absorbed in H2O)*P(not absorbed in acrylic)
@@ -313,7 +324,7 @@ static double get_expected_charge(double x, double beta, double theta0, double *
 
     *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
 
-    return (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
+    *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
 }
 
 double time_cdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma)
@@ -403,50 +414,48 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
 static void integrate_path_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma)
 {
     size_t i;
-    static double qs[MAX_NPOINTS];
-    static double q_indirects[MAX_NPOINTS];
-    static double ts[MAX_NPOINTS];
-    static double ts2[MAX_NPOINTS];
-    double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o, q_delta_ray, q_indirect_delta_ray, dx, tmp;
-
-    if (n > MAX_NPOINTS) {
-        fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n");
-        exit(1);
-    }
-
-    /* FIXME: I just calculate delta assuming 400 nm light. */
-    wavelength0 = 400.0;
-    n_d2o = get_index_snoman_d2o(wavelength0);
-    n_h2o = get_index_snoman_h2o(wavelength0);
+    double pos[3], t, q, r, qd, rd, dx, tmp, q_sum, r_sum, t_sum, t2_sum;
 
+    q_sum = 0.0;
+    r_sum = 0.0;
+    t_sum = 0.0;
+    t2_sum = 0.0;
     for (i = 0; i < n; i++) {
         pos[0] = pos0[0] + x[i]*dir0[0];
         pos[1] = pos0[1] + x[i]*dir0[1];
         pos[2] = pos0[2] + x[i]*dir0[2];
 
-        get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);
+        get_expected_charge_shower(p, pos, dir0, pmt, &q, &r, &qd, &rd, &t);
 
-        t = x[i]/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
+        q = q*pdf[i] + qd/p->range;
+        r = r*pdf[i] + rd/p->range;
+        t += x[i]/SPEED_OF_LIGHT;
 
-        qs[i] = get_expected_charge_shower(p, pos, dir0, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o, &q_delta_ray, &q_indirect_delta_ray);
-        qs[i] = qs[i]*pdf[i] + q_delta_ray/p->range;
-        q_indirects[i] = q_indirects[i]*pdf[i] + q_indirect_delta_ray/p->range;
-        ts[i] = t*qs[i];
-        ts2[i] = t*t*qs[i];
+        if (i == 0 || i == (n - 1)) {
+            q_sum += q;
+            r_sum += r;
+            t_sum += t*q;
+            t2_sum += t*t*q;
+        } else {
+            q_sum += 2*q;
+            r_sum += 2*r;
+            t_sum += 2*t*q;
+            t2_sum += 2*t*t*q;
+        }
     }
 
     dx = x[1]-x[0];
-    *mu_direct = trapz(qs,dx,n);
-    *mu_indirect = trapz(q_indirects,dx,n);
+    *mu_direct = q_sum*dx*0.5;
+    *mu_indirect = r_sum*dx*0.5;
 
     if (*mu_direct == 0.0) {
         *time = 0.0;
         *sigma = PMT_TTS;
     } else {
-        *time = trapz(ts,dx,n)/(*mu_direct);
+        *time = t_sum*dx*0.5/(*mu_direct);
 
         /* Variance in the time = E(t^2) - E(t)^2. */
-        tmp = trapz(ts2,dx,n)/(*mu_direct) - (*time)*(*time);
+        tmp = t2_sum*dx*0.5/(*mu_direct) - (*time)*(*time);
 
         if (tmp >= 0) {
             *sigma = sqrt(tmp + PMT_TTS*PMT_TTS);
@@ -461,24 +470,12 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0
 static void integrate_path(path *p, int pmt, double *mu_direct, double *mu_indirect, double *time)
 {
     size_t i;
-    static double qs[MAX_NPOINTS];
-    static double q_indirects[MAX_NPOINTS];
-    static double ts[MAX_NPOINTS];
-    double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o, x;
-
-    if (p->len > MAX_NPOINTS) {
-        fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n");
-        exit(1);
-    }
-
-    /* FIXME: I just calculate delta assuming 400 nm light. */
-    wavelength0 = 400.0;
-    n_d2o = get_index_snoman_d2o(wavelength0);
-    n_h2o = get_index_snoman_h2o(wavelength0);
+    double dir[3], pos[3], t, theta0, beta, q, r, q_sum, r_sum, t_sum, dx;
 
+    q_sum = 0.0;
+    r_sum = 0.0;
+    t_sum = 0.0;
     for (i = 0; i < p->len; i++) {
-        x = p->s[i];
-
         pos[0] = p->x[i];
         pos[1] = p->y[i];
         pos[2] = p->z[i];
@@ -489,30 +486,35 @@ static void integrate_path(path *p, int pmt, double *mu_direct, double *mu_indir
 
         beta = p->beta[i];
 
-        t = p->t[i];
-
         if (p->n > 0) {
             theta0 = p->theta0;
         } else {
-            theta0 = fmax(MIN_THETA0,p->theta0*sqrt(x));
+            theta0 = fmax(MIN_THETA0,p->theta0*sqrt(p->s[i]));
             theta0 = fmin(MAX_THETA0,theta0);
         }
 
-        get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);
+        get_expected_charge(beta, theta0, pos, dir, pmt, &q, &r, &t);
 
-        t += l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
+        t += p->t[i];
 
-        qs[i] = get_expected_charge(x, beta, theta0, pos, dir, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o);
-        ts[i] = t*qs[i];
+        if (i == 0 || i == p->len - 1) {
+            q_sum += q;
+            r_sum += r;
+            t_sum += t*q;
+        } else {
+            q_sum += 2*q;
+            r_sum += 2*r;
+            t_sum += 2*t*q;
+        }
     }
 
-    double dx = p->s[1]-p->s[0];
-    *mu_direct = trapz(qs,dx,p->len);
-    *mu_indirect = trapz(q_indirects,dx,p->len);
-    *time = trapz(ts,dx,p->len);
+    dx = p->s[1] - p->s[0];
+    *mu_direct = q_sum*dx*0.5;
+    *mu_indirect = r_sum*dx*0.5;
+    *time = t_sum*dx*0.5;
 }
 
-static double get_total_charge_approx(double beta0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected, double n_d2o, double n_h2o, double cos_theta_cerenkov, double sin_theta_cerenkov)
+static double get_total_charge_approx(double beta0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected, double cos_theta_cerenkov, double sin_theta_cerenkov)
 {
     /* Returns the approximate expected number of photons seen by PMT `i` using
      * an analytic formula.
@@ -607,7 +609,7 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
 
     *t = 0.0;
     *mu_reflected = 0.0;
-    if (beta < 1/n_d2o) return 0.0;
+    if (beta < 1/avg_index_d2o) return 0.0;
 
     /* `prob` is the number of photons emitted per cm by the particle at a
      * distance `s` along the track. */
@@ -638,7 +640,7 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
 
     theta0 = fmax(theta0*sqrt(s),MIN_THETA0);
 
-    frac = sqrt(2)*n_d2o*x*beta0*theta0;
+    frac = sqrt(2)*avg_index_d2o*x*beta0*theta0;
 
     f = get_weighted_pmt_response(-cos_theta_pmt);
     f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
@@ -649,15 +651,15 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
     get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o);
 
     /* Assume the particle is travelling at the speed of light. */
-    *t = s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
+    *t = s/SPEED_OF_LIGHT + l_d2o*avg_index_d2o/SPEED_OF_LIGHT + l_h2o*avg_index_h2o/SPEED_OF_LIGHT;
 
-    charge = n_d2o*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n_d2o*(smax-z)*beta0)/frac) + erf((-a+b*s+n_d2o*z*beta0)/frac))/(b+n_d2o*beta0)/(4*M_PI);
+    charge = avg_index_d2o*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+avg_index_d2o*(smax-z)*beta0)/frac) + erf((-a+b*s+avg_index_d2o*z*beta0)/frac))/(b+avg_index_d2o*beta0)/(4*M_PI);
 
     /* Add expected number of photons from electromagnetic shower. */
     if (p->shower_photons > 0)
-        charge += get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/n_d2o)*omega/(2*M_PI);
+        charge += get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/avg_index_d2o)*omega/(2*M_PI);
     if (p->delta_ray_photons > 0)
-        charge += get_weighted_quantum_efficiency()*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1.0/n_d2o)*omega/(2*M_PI);
+        charge += get_weighted_quantum_efficiency()*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1.0/avg_index_d2o)*omega/(2*M_PI);
 
     *mu_reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
 
@@ -729,7 +731,7 @@ static int get_smax(particle *p, double beta_min, double range, double *smax)
     return status;
 }
 
-static double guess_time(double *pos, double *dir, double *pmt_pos, double smax, double n_d2o, double n_h2o, double cos_theta_cerenkov, double sin_theta_cerenkov)
+static double guess_time(double *pos, double *dir, double *pmt_pos, double smax, double cos_theta_cerenkov, double sin_theta_cerenkov)
 {
     /* Returns an approximate time at which a PMT is most likely to get hit
      * from Cerenkov light.
@@ -777,7 +779,7 @@ static double guess_time(double *pos, double *dir, double *pmt_pos, double smax,
     get_path_length(tmp,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);
 
     /* Assume the particle is travelling at the speed of light. */
-    return s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
+    return s/SPEED_OF_LIGHT + l_d2o*avg_index_d2o/SPEED_OF_LIGHT + l_h2o*avg_index_h2o/SPEED_OF_LIGHT;
 }
 
 static double getKineticEnergy(double x, void *p)
@@ -899,7 +901,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
     static double logp[MAX_PE], nll[MAX_PMTS];
     double range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
     particle *p;
-    double wavelength0, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov;
+    double cos_theta_cerenkov, sin_theta_cerenkov;
     double logp_path;
     double mu_reflected;
     path *path;
@@ -1017,10 +1019,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
          *
          * These are computed at the beginning of this function and then passed to
          * the different functions to avoid recomputing them on the fly. */
-        wavelength0 = 400.0;
-        n_d2o = get_index_snoman_d2o(wavelength0);
-        n_h2o = get_index_snoman_h2o(wavelength0);
-        cos_theta_cerenkov = 1/(n_d2o*beta0);
+        cos_theta_cerenkov = 1/(avg_index_d2o*beta0);
         sin_theta_cerenkov = sqrt(1-pow(cos_theta_cerenkov,2));
 
         mu_indirect[j] = 0.0;
@@ -1032,7 +1031,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
             if (hit_only && !ev->pmt_hits[i].hit) continue;
 
             if (fast) {
-                mu_direct[i][j] = get_total_charge_approx(beta0, v[j].pos, v[j].dir, p, i, smax, theta0, &ts[i][j], &mu_reflected, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov);
+                mu_direct[i][j] = get_total_charge_approx(beta0, v[j].pos, v[j].dir, p, i, smax, theta0, &ts[i][j], &mu_reflected, cos_theta_cerenkov, sin_theta_cerenkov);
                 ts[i][j] += v[j].t0;
                 mu_indirect[j] += mu_reflected;
                 continue;
@@ -1055,7 +1054,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
                  * I should test it to see if we can get away with something even
                  * simpler (like just computing the distance from the start of the
                  * track to the PMT). */
-                ts[i][j] = v[j].t0 + guess_time(v[j].pos,v[j].dir,pmts[i].pos,smax,n_d2o,n_h2o,cos_theta_cerenkov,sin_theta_cerenkov);
+                ts[i][j] = v[j].t0 + guess_time(v[j].pos,v[j].dir,pmts[i].pos,smax,cos_theta_cerenkov,sin_theta_cerenkov);
             }
 
             integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu_shower[i][j],&q_indirect,&result,&ts_sigma[i][j]);
diff --git a/src/optics.c b/src/optics.c
index 830ada9..6b6d761 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -29,6 +29,8 @@
 
 char optics_err[256];
 
+double avg_index_h2o, avg_index_d2o;
+
 /* Absorption coefficients for H2O, D2O, and acrylic as a function of
  * wavelength from SNOMAN.
  *
@@ -425,6 +427,28 @@ static double gsl_absorption_length_snoman_acrylic(double wavelength, void *para
     return qe*exp(-x/get_absorption_length_snoman_acrylic(wavelength))/pow(wavelength,2);
 }
 
+/* Returns the product of the index of refraction in D2O multiplied by the
+ * quantum efficiency and the Cerenkov spectrum. */
+static double gsl_index_d2o(double wavelength, void *params)
+{
+    double qe;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe*get_index_snoman_d2o(wavelength)/pow(wavelength,2);
+}
+
+/* Returns the product of the index of refraction in H2O multiplied by the
+ * quantum efficiency and the Cerenkov spectrum. */
+static double gsl_index_h2o(double wavelength, void *params)
+{
+    double qe;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe*get_index_snoman_h2o(wavelength)/pow(wavelength,2);
+}
+
 static double gsl_cerenkov(double wavelength, void *params)
 {
     /* Returns the quantum efficiency multiplied by the Cerenkov spectrum. */
@@ -458,6 +482,30 @@ int optics_init(dict *db)
         return -1;
     }
 
+    F.function = &gsl_index_d2o;
+    F.params = NULL;
+
+    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+    if (status) {
+        sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+        return -1;
+    }
+
+    avg_index_d2o = result/norm;
+
+    F.function = &gsl_index_h2o;
+    F.params = NULL;
+
+    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+    if (status) {
+        sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+        return -1;
+    }
+
+    avg_index_h2o = result/norm;
+
     for (i = 0; i < LEN(fabs_x); i++) {
         fabs_x[i] = xlo + (xhi-xlo)*i/(LEN(fabs_x)-1);
     }
diff --git a/src/optics.h b/src/optics.h
index 70ff53b..4b41c6b 100644
--- a/src/optics.h
+++ b/src/optics.h
@@ -22,6 +22,8 @@
 /* Global error string when optics_init() returns -1. */
 extern char optics_err[256];
 
+extern double avg_index_h2o, avg_index_d2o;
+
 /* Initialize the optics data by reading in the RSPR bank and precomputing the
  * average absorption and scattering tables. */
 int optics_init(dict *db);