--- ray/src/px/pf3.c	1992/09/21 12:14:12	2.2
+++ ray/src/px/pf3.c	1993/06/18 13:36:20	2.5
@@ -12,6 +12,8 @@ static char SCCSid[] = "$SunId$ LBL";
 
 #include  <stdio.h>
 
+#include  <math.h>
+
 #include  "color.h"
 
 #define	 FTINY		1e-7
@@ -36,24 +38,27 @@ extern COLOR  *scanout;		/* output scan line */
 
 extern char  *progname;
 
-float  *exptable;		/* exponent table */
+float  *gausstable;		/* gauss lookup table */
 
-#define	 lookexp(x)		exptable[(int)(-10.*(x)+.5)]
+#define	 lookgauss(x)		gausstable[(int)(-10.*(x)+.5)]
 
 
 initmask()			/* initialize gaussian lookup table */
 {
 	extern char  *malloc();
-	extern double  exp();
+	double  d;
 	register int  x;
 
-	exptable = (float *)malloc(100*sizeof(float));
-	if (exptable == NULL) {
+	gausstable = (float *)malloc(100*sizeof(float));
+	if (gausstable == NULL) {
 		fprintf(stderr, "%s: out of memory in initmask\n", progname);
 		quit(1);
 	}
-	for (x = 0; x < 100; x++)
-		exptable[x] = exp(-x*0.1);
+	d = x_c*y_r*0.25/rad/rad;
+	gausstable[0] = exp(-d)/sqrt(d);
+	for (x = 1; x < 100; x++)
+		if ((gausstable[x] = exp(-x*0.1)/sqrt(x*0.1)) > gausstable[0])
+			gausstable[x] = gausstable[0];
 }
 
 
@@ -71,16 +76,16 @@ int  c, r;
 	wsum = 0;
 	setcolor(csum, 0.0, 0.0, 0.0);
 	for (y = ycent+1-yrad; y <= ycent+yrad; y++) {
-		if (y < 0 || y >= yres)
-			continue;
-		d = y_r < 1.0 ? y_r*y - r : y - ycent;
+		if (y < 0) continue;
+		if (y >= yres) break;
+		d = y_r < 1.0 ? y_r*y - r : (double)(y - ycent);
 		if (d > 0.5+FTINY || d < -0.5-FTINY)
 			continue;
 		scan = scanin[y%barsize];
 		for (x = xcent+1-xrad; x <= xcent+xrad; x++) {
-			if (x < 0 || x >= xres)
-				continue;
-			d = x_c < 1.0 ? x_c*x - c : x - xcent;
+			if (x < 0) continue;
+			if (x >= xres) break;
+			d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent);
 			if (d > 0.5+FTINY || d < -0.5-FTINY)
 				continue;
 			wsum++;
@@ -106,15 +111,15 @@ int  c, r;
 	wsum = FTINY;
 	setcolor(csum, 0.0, 0.0, 0.0);
 	for (y = ycent-yrad; y <= ycent+yrad; y++) {
-		if (y < 0 || y >= yres)
-			continue;
+		if (y < 0) continue;
+		if (y >= yres) break;
 		dy = (y_r*(y+.5) - (r+.5))/rad;
 		scan = scanin[y%barsize];
 		for (x = xcent-xrad; x <= xcent+xrad; x++) {
-			if (x < 0 || x >= xres)
-				continue;
+			if (x < 0) continue;
+			if (x >= xres) break;
 			dx = (x_c*(x+.5) - (c+.5))/rad;
-			weight = lookexp(-(dx*dx + dy*dy));
+			weight = lookgauss(-(dx*dx + dy*dy));
 			wsum += weight;
 			copycolor(ctmp, scan[x]);
 			scalecolor(ctmp, weight);