diff --git a/jdmath/src/dinterpo.c b/jdmath/src/dinterpo.c
index b1d9765..15e5af7 100644
--- a/jdmath/src/dinterpo.c
+++ b/jdmath/src/dinterpo.c
@@ -44,32 +44,38 @@ unsigned int JDMbinary_search_d (double x, double *xp, unsigned int n)
 	  }
 	else n0 = n2;
      }
+   if (x >= xp[n0]) return n1;
    return n0;
 }
 
 double JDMinterpolate_d (double x, double *xp, double *yp, unsigned int n)
 {
-   unsigned int n0, n1;
-   double x0, x1;
-
-   n0 = JDMbinary_search_d (x, xp, n);
-
-   x0 = xp[n0];
-   n1 = n0 + 1;
-
-   if (x == x0)
-     return yp[n0];
-   if (n1 == n)
-     {
-	if (n0 == 0)
-	  return yp[n0];
-	n1 = n0 - 1;
-     }
-
-   x1 = xp[n1];
-   if (x1 == x0) return yp[n0];
-
-   return yp[n0] + (yp[n1] - yp[n0]) / (x1 - x0) * (x - x0);
+	unsigned int n0, n1;
+	double x0, x1;
+
+	if (n == 1)
+		return yp[0];
+	n1 = JDMbinary_search_f(x, xp, n);
+	n0 = n1 - 1;
+
+	if (x == xp[n1])
+		return yp[n1];
+	if (n1 == n)
+	{
+		n1--;
+		n0--;
+	}
+	if (n1 == 0)
+	{
+		n0 = 1;
+	}
+
+	x0 = xp[n0];
+	x1 = xp[n1];
+	if (x1 == x0)
+		return yp[n1];
+
+	return yp[n0] + (yp[n1] - yp[n0]) / (x1 - x0) * (x - x0);
 }
 
 int JDMinterpolate_dvector (double *xp, double *yp, unsigned int n,
diff --git a/jdmath/src/finterpo.c b/jdmath/src/finterpo.c
index 74aaf6e..ec84f80 100644
--- a/jdmath/src/finterpo.c
+++ b/jdmath/src/finterpo.c
@@ -27,6 +27,13 @@
 
 #include "jdmath.h"
 
+/*
+ * Binary Search
+ *
+ * This routine returns the index of the first element in the array xp
+ * that is greater than or equal to x.  If all elements are less than x,
+ * then n is returned.
+ */
 unsigned int JDMbinary_search_f (float x, float *xp, unsigned int n)
 {
    unsigned int n0, n1, n2;
@@ -42,10 +49,7 @@ unsigned int JDMbinary_search_f (float x, float *xp, unsigned int n)
 	     if (xp[n2] == x) return n2;
 	     n1 = n2;
 	  }
-	else
-	  {
-	     n0 = n2;
-	  }
+	else n0 = n2;
      }
 
    if (x >= xp[n0]) return n1;
@@ -57,22 +61,25 @@ float JDMinterpolate_f (float x, float *xp, float *yp, unsigned int n)
    unsigned int n0, n1;
    double x0, x1;
 
-   n0 = JDMbinary_search_f (x, xp, n);
-
-   x0 = xp[n0];
-   n1 = n0 + 1;
+   if (n == 1) return yp[0];
+   n1 = JDMbinary_search_f (x, xp, n);
+   n0 = n1 - 1;
 
-   if (x == x0)
-     return yp[n0];
+   if (x == xp[n1])
+	   return yp[n1];
    if (n1 == n)
      {
-	if (n0 == 0)
-	  return yp[n0];
-	n1 = n0 - 1;
+	    n1--;
+		n0--;
      }
+   if (n1 == 0)
+	 {
+		 n0 = 1;
+	 }
 
+   x0 = xp[n0];
    x1 = xp[n1];
-   if (x1 == x0) return yp[n0];
+   if (x1 == x0) return yp[n1];
 
    return yp[n0] + (yp[n1] - yp[n0]) / (x1 - x0) * (x - x0);
 }
@@ -192,14 +199,15 @@ float JDMlog_interpolate_f (float x, float *xp, float *yp, unsigned int n)
    unsigned int n0, n1;
    double x0, x1;
 
-   n0 = JDMbinary_search_f (x, xp, n);
+   n1 = JDMbinary_search_f (x, xp, n);
+   n0 = n1 - 1;
 
-   x0 = xp[n0];
-   n1 = n0 + 1;
-
-   if ((x == x0) || (n1 == n)) return yp[n0];
+   if (n1 == 0) return yp[0];
+   if (n1 == n) return yp[n - 1];
 
+   x0 = xp[n0];
    x1 = xp[n1];
+   if (x == x1) return yp[n1];
    if (x1 == x0) return yp[n0];
 
    return yp[n0] + (yp[n1] - yp[n0]) * (log(x/x0) / log (x1/x0));
diff --git a/jdmath/src/hist.c b/jdmath/src/hist.c
index ab72ce8..1afb58a 100644
--- a/jdmath/src/hist.c
+++ b/jdmath/src/hist.c
@@ -40,27 +40,6 @@
 #define FLOAT_TYPE float
 #include "hist.c"
 #else
-#if 0
-static unsigned int binary_search_d (double x, double *xp, unsigned int n)
-{
-   unsigned int n0, n1, n2;
-
-   n0 = 0;
-   n1 = n;
-
-   while (n1 > n0 + 1)
-     {
-	n2 = (n0 + n1) / 2;
-	if (xp[n2] >= x)
-	  {
-	     if (xp[n2] == x) return n2;
-	     n1 = n2;
-	  }
-	else n0 = n2;
-     }
-   return n0;
-}
-#endif
 /*
  * histogram routines
  *
@@ -68,6 +47,7 @@ static unsigned int binary_search_d (double x, double *xp, unsigned int n)
  * of values y_i to be grouped into the histogram.  The bin size of the
  * nth bin is given by x_{i+1} - x_i, except for the last bin, which is
  * assumed to be of infinite width.
+ * This means that values of pts that are less than x_0 are ignored.
  */
 
 /* If reverse_indices is NON-NULL, it is assumed to point to an array of
@@ -102,10 +82,10 @@ int JDMHISTOGRAM_FUNCTION (FLOAT_TYPE *pts, unsigned int npts,
 	  continue;
 
 	j = BINARY_SEARCH (val, bin_edges, nbins);
-	histogram[j] += 1;
+	histogram[j - 1] += 1;
 
 	if (reverse_indices != NULL)
-	  reverse_indices[i] = (int) j;
+	  reverse_indices[i] = (int) j - 1;
      }
    return 0;
 }
diff --git a/marx/doc/examples/user-source/eventlist.c b/marx/doc/examples/user-source/eventlist.c
index 831f40d..d6c9c99 100644
--- a/marx/doc/examples/user-source/eventlist.c
+++ b/marx/doc/examples/user-source/eventlist.c
@@ -594,7 +594,9 @@ static double map_pi_to_line (short pi)
 	return b->energies[0];
      }
    i = JDMbinary_search_f (JDMrandom (), b->cum_strengths, b->num_energies);
-   return b->energies[i];
+   if (i == 0)
+     return b->energies[0];
+   return b->energies[i - 1];
 }
 
 static int read_data (double *tp, double *xp, double *yp, short *chanp)
diff --git a/marx/libsrc/acis_fef.c b/marx/libsrc/acis_fef.c
index bf255e2..cecbfec 100644
--- a/marx/libsrc/acis_fef.c
+++ b/marx/libsrc/acis_fef.c
@@ -553,19 +553,20 @@ static int normalize_gaussians (Gauss_Parm_Type *gaussians, unsigned int num,
 	/* return 0; */
      }
 
-   if (total_neg_area != 0.0)
-     flags |= HAS_NEG_AMP_GAUSSIANS;
-
-   g = gaussians;
-   if (total_pos_area > 0) while (g < gmax)
-     {
-	/* Since we interpolate over pairs of gaussians, do not normalize
-	 * each member of the pair separately.
-	 */
-	/* g->amp /= total_pos_area; */
-	g->cum_area /= total_pos_area;
-	g++;
-     }
+     if (total_neg_area != 0.0)
+       flags |= HAS_NEG_AMP_GAUSSIANS;
+
+     g = gaussians;
+     if (total_pos_area > 0)
+       while (g < gmax)
+       {
+         /* Since we interpolate over pairs of gaussians, do not normalize
+          * each member of the pair separately.
+          */
+         /* g->amp /= total_pos_area; */
+         g->cum_area /= total_pos_area;
+         g++;
+       }
 
    total_area = total_pos_area - total_neg_area;
    if (total_area != 0)
@@ -990,18 +991,22 @@ int marx_apply_acis_rmf (int ccd_id, float x, float y,
      }
 
    num_energies = f->num_energies;
-   i = JDMbinary_search_f (energy, f->energies, num_energies);
-   j = i + 1;
+   j = JDMbinary_search_f (energy, f->energies, num_energies);
+   if (j == 0)
+   {
+     /* extrapolate */
+     j++;
+   }
    if (j == num_energies)
-     {
-	/* extrapolate */
-	i--;
-	j--;
-     }
+   {
+	   /* extrapolate */
+     j--;
+   }
 
-   t = (energy - f->energies[i])/(f->energies[j] - f->energies[i]);
+   t = (energy - f->energies[j-1])/(f->energies[j] - f->energies[j-1]);
    g0 = f->gaussians + i * num_gaussians;
    g1 = g0 + num_gaussians;
+   //marx_message("acis_fef: t = %f, energy= %f, f->energies[j-1] = %f, f->energies[j] = %f \n", t, energy, f->energies[j-1], f->energies[j]);
 
    for (i = 0; i < num_gaussians; i++)
      {
@@ -1031,10 +1036,10 @@ int marx_apply_acis_rmf (int ccd_id, float x, float y,
 
    if (flags & HAS_TOTAL_NEG_AREA)
      {
-	marx_message ("Warning occurred for ccdid=%d,chipx=%f,chipy=%f,energy=%f\n",
-		      ccd_id, x, y, energy);
-	return -1;
-	/* flags &= ~HAS_TOTAL_NEG_AREA; */
+       marx_message("Warning occurred for ccdid=%d,chipx=%f,chipy=%f,energy=%f,flags=%d,flagsandneg=%d\n",
+                    ccd_id, x, y, energy, flags, flags & HAS_TOTAL_NEG_AREA);
+       return -1;
+       /* flags &= ~HAS_TOTAL_NEG_AREA; */
      }
 
    if (flags == 0)
diff --git a/marx/libsrc/acis_subpix.c b/marx/libsrc/acis_subpix.c
index a28d312..8f182a8 100644
--- a/marx/libsrc/acis_subpix.c
+++ b/marx/libsrc/acis_subpix.c
@@ -284,7 +284,7 @@ int marx_compute_acis_subpix (Marx_Subpix_Table_Type *stt,
 			      int ccd, float energy, int fltgrade, float *dxp, float *dyp)
 {
    Subpix_Type *s;
-   unsigned int i, j, n;
+   unsigned int j, n;
    double w0, w1;
 
    if (stt == NULL)
@@ -301,18 +301,22 @@ int marx_compute_acis_subpix (Marx_Subpix_Table_Type *stt,
      }
 
    n = (unsigned int) s->num_energies;
-   i = JDMbinary_search_f (energy, s->energies, n);
-   j = i+1;
+   j = JDMbinary_search_f (energy, s->energies, n);
+   if (j == 0)
+   {
+     /* extrapolate */
+     j++;
+   }
    if (j == n)
-     {
-	j = i;
-	i--;
-     }
+   {
+     /* extrapolate */
+     j--;
+   }
    // Linear interpolation on the energy grid to get dxp, dyp
-   w1 = ((double)energy - s->energies[i])/(s->energies[j] - s->energies[i]);
+   w1 = ((double)energy - s->energies[j - 1])/(s->energies[j] - s->energies[j - 1]);
    w0 = 1.0 - w1;
-   *dxp = w0 * s->dxs[i] + w1 * s->dxs[j];
-   *dyp = w0 * s->dys[i] + w1 * s->dys[j];
+   *dxp = w0 * s->dxs[j - 1] + w1 * s->dxs[j];
+   *dyp = w0 * s->dys[j - 1] + w1 * s->dys[j];
 
    return 0;
 }
diff --git a/marx/libsrc/aciscontam.c b/marx/libsrc/aciscontam.c
index 14b7d73..e6561e5 100644
--- a/marx/libsrc/aciscontam.c
+++ b/marx/libsrc/aciscontam.c
@@ -240,18 +240,26 @@ static void compute_tau_values (Single_Component_Contam_Type *contam,
    double t = _Marx_TStart_MJDsecs;
    double t0, t1;
 
-   n0 = JDMbinary_search_d (t, times, ntimes);
-   n1 = n0+1;
+   n1 = JDMbinary_search_d (t, times, ntimes);
+   n0 = n1 - 1;
 
-   if (n1 == ntimes)
-     {
-	if (n0 == 0)
+   if (ntimes == 1)
 	  {
 	     contam->tau_0s[layer] = tau0s[0];
 	     contam->tau_1s[layer] = tau1s[0];
 	     return;
 	  }
-	n1 = n0-1;
+
+  /* Out of range. Extrapolate. */
+  if (n1 == ntimes)
+     {
+      n1--;
+      n0--;
+     }
+  if (n1 == 0)
+     {
+      n1 = 1;
+      n0 = 0;
      }
 
    t0 = times[n0];
diff --git a/marx/libsrc/diffract.c b/marx/libsrc/diffract.c
index 91dd18b..70847ff 100644
--- a/marx/libsrc/diffract.c
+++ b/marx/libsrc/diffract.c
@@ -337,9 +337,15 @@ double marx_compute_grating_efficiency (double energy, int order,
      }
 
    /* Find where energy lies on the Energies grid. */
-   n0 = JDMbinary_search_f (energy, Energies, Num_Energies);
-   n1 = n0 + 1;
-
+   n1 = JDMbinary_search_f (energy, Energies, Num_Energies);
+   n0 = n1 - 1;
+   if (n1 == 0)
+  {
+    /* Oops.  energy is out of range.
+     * Use last first points to extrapolate. */
+    n0++;
+    n1++;
+  }
    if (n1 == Num_Energies)
      {
 	/* Oops.  energy is out of range.
@@ -356,7 +362,7 @@ double marx_compute_grating_efficiency (double energy, int order,
 
    e1 = Energies[n1];
    eff1 = compute_efficiency (order, n1, geom);
-
+   marx_message ("energy = %g, e0 = %g, e1 = %g\n", energy, e0, e1);
    return eff0 + (eff1 - eff0)/(e1 - e0) * (energy - e0);
 }
 
@@ -733,6 +739,14 @@ static int diffract_photon (Grating_Type *g, double theta,
    d = JDMv_cross_prod (n, l);
 
    if (gs != NULL)
+     /* HMG 11/17/2023:
+      * As far as I can tell, this code is not used at the moment.
+      * It's implemented here so that we can switch it on with detailed
+      * sector information, however, that information has never been obtained
+      * in calibration and thus is not available in MARX/CALDB.
+      * So, gs (the variable holding the grating sector information) is always NULL
+      * and the first branch of the if is never used.
+      */
      {
 	unsigned int sector_num;
 	double sector;
@@ -742,7 +756,8 @@ static int diffract_photon (Grating_Type *g, double theta,
 	if (sector < 0) sector = 2*PI + sector;
 
 	sector_num = JDMbinary_search_d (sector, gs->min_angle, gs->num_sectors);
-	if ((gs->max_angle[sector_num] <= sector)
+  sector_num = sector_num - 1;
+  if ((gs->max_angle[sector_num] <= sector)
 	    || (gs->min_angle[sector_num] > sector))
 	  return -1;
 
diff --git a/marx/libsrc/mblur.c b/marx/libsrc/mblur.c
index 0521daa..b09bb17 100644
--- a/marx/libsrc/mblur.c
+++ b/marx/libsrc/mblur.c
@@ -234,10 +234,20 @@ int _marx_mirror_blur (Marx_Photon_Type *pt) /*{{{*/
 	at = photon_attributes + sorted_index[i];
 
 	energy = (float) at->energy;
-	/* Now find which two energy curves to use: j and j + 1 */
+	/* Now find which two energy curves to use: j - 1 and j */
 	j = JDMbinary_search_f (energy, energies, num_energies);
-	e0 = energies[j];
-	e1 = energies[j + 1];
+  if (j == 0)
+  {
+    /* extrapolate */
+    j++;
+  }
+  if (j == num_energies)
+  {
+    /* extrapolate */
+    j--;
+  }
+  e0 = energies[j - 1];
+	e1 = energies[j];
 
 	/* Ideally a random number between 0 and 1 would be used.  However,
 	 * the blur data close to 1.0 cannot be trusted.
@@ -246,10 +256,10 @@ int _marx_mirror_blur (Marx_Photon_Type *pt) /*{{{*/
 
 	/* Now interpolate thetas on the two energy curves */
 	theta0 = JDMinterpolate_f (r, encircled_energies,
-				   Mirr_Blur.thetas[at->mirror_shell][j],
+				   Mirr_Blur.thetas[at->mirror_shell][j - 1],
 				   num_ee);
 	theta1 = JDMinterpolate_f (r, encircled_energies,
-				   Mirr_Blur.thetas[at->mirror_shell][j + 1],
+				   Mirr_Blur.thetas[at->mirror_shell][j],
 				   num_ee);
 	blur = theta0 + (theta1 - theta0) * (energy - e0) / (e1 - e0);
 
diff --git a/marx/libsrc/wfold.c b/marx/libsrc/wfold.c
index 423420e..fb325b8 100644
--- a/marx/libsrc/wfold.c
+++ b/marx/libsrc/wfold.c
@@ -347,16 +347,21 @@ double marx_wfold_table_interp (Marx_WFold_Table_Type *table, double energy, dou
    e_alpha = energy * sin_alpha;
    i = JDMbinary_search_d (e_alpha, table->e_alphas, table->num_arrays);
 
-   if (i + 1 == table->num_arrays)
+   if (i == table->num_arrays)
      {
-	/* Back off one and extrapolate outside range */
-	i--;
+	      /* Back off one and extrapolate outside range */
+	      i--;
+     }
+     if (i == 0)
+     {
+       /* Back off one and extrapolate outside range */
+       i++;
      }
 
-   theta_0 = interpolate_theta (table->fold_arrays[i], r);
-   theta_1 = interpolate_theta (table->fold_arrays[i + 1], r);
-   e_alpha_0 = table->e_alphas[i];
-   e_alpha_1 = table->e_alphas[i + 1];
+   theta_0 = interpolate_theta (table->fold_arrays[i - 1], r);
+   theta_1 = interpolate_theta (table->fold_arrays[i], r);
+   e_alpha_0 = table->e_alphas[i - 1];
+   e_alpha_1 = table->e_alphas[i];
    theta = theta_0 + (theta_1 - theta_0) * (e_alpha - e_alpha_0) / (e_alpha_1 - e_alpha_0);
 
    /* Arc seconds??  Assume radians */