aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/likelihood.c181
-rw-r--r--src/optics.c48
-rw-r--r--src/optics.h2
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);
generated by cgit (git 2.34.1) at 2025-05-06 21:35:23 +0000