--- ray/src/common/interp2d.c	2013/02/11 23:33:35	2.5
+++ ray/src/common/interp2d.c	2013/02/12 17:47:58	2.7
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: interp2d.c,v 2.5 2013/02/11 23:33:35 greg Exp $";
+static const char RCSid[] = "$Id: interp2d.c,v 2.7 2013/02/12 17:47:58 greg Exp $";
 #endif
 /*
  * General interpolation method for unstructured values on 2-D plane.
@@ -17,7 +17,7 @@ static const char RCSid[] = "$Id: interp2d.c,v 2.5 201
  * each of NI2DIR surrounding directions.  To speed this
  * calculation, we sort the data into half-planes and apply
  * simple tests to see which neighbor is closest in each
- * direction.  Once we have our approximate neighborhood
+ * angular slice.  Once we have our approximate neighborhood
  * for a sample, we can use it in a modified Gaussian weighting
  * with allowing local anisotropy.  Harmonic weighting is added
  * to reduce the influence of distant neighbors.  This yields a
@@ -90,7 +90,7 @@ interp2_spacing(INTERP2 *ip, double mind)
 {
 	if (mind <= 0)
 		return;
-	if ((.998*ip->dmin <= mind) && (mind <= 1.002*ip->dmin))
+	if ((.998*ip->dmin <= mind) & (mind <= 1.002*ip->dmin))
 		return;
 	if (ip->da != NULL) {	/* will need to recompute distribution */
 		free(ip->da);
@@ -238,7 +238,7 @@ interp2_analyze(INTERP2 *ip)
 static double
 get_wt(const INTERP2 *ip, const int i, double x, double y)
 {
-	double	dir, rd, d2;
+	double	dir, rd, r2, d2;
 	int	ri;
 				/* get relative direction */
 	x -= ip->spt[i][0];
@@ -251,9 +251,12 @@ get_wt(const INTERP2 *ip, const int i, double x, doubl
 	rd -= (double)ri;
 	rd = (1.-rd)*ip->da[i][ri] + rd*ip->da[i][(ri+1)%NI2DIR];
 	rd = ip->smf * DECODE_DIA(ip, rd);
+	r2 = 2.*rd*rd;
 	d2 = x*x + y*y;
+	if (d2 > 21.*r2)	/* result would be < 1e-9 */
+		return(.0);
 				/* Gaussian times harmonic weighting */
-	return( exp(d2/(-2.*rd*rd)) * ip->dmin/(ip->dmin + sqrt(d2)) );
+	return( exp(-d2/r2) * ip->dmin/(ip->dmin + sqrt(d2)) );
 }
 
 /* Assign full set of normalized weights to interpolate the given position */