--- ray/src/rt/ambcomp.c	1991/10/22 12:15:41	1.15
+++ ray/src/rt/ambcomp.c	1996/04/23 11:06:54	2.6
@@ -27,9 +27,7 @@ typedef struct {
 	short  nt, np;		/* number of theta and phi directions */
 }  AMBHEMI;		/* ambient sample hemisphere */
 
-extern double  sin(), cos(), sqrt();
 
-
 static int
 ambcmp(d1, d2)				/* decreasing order */
 AMBSAMP  *d1, *d2;
@@ -141,14 +139,15 @@ FVECT  pg, dg;
 			dp->n = 0;
 			if (divsample(dp, &hemi, r) < 0)
 				goto oopsy;
+			arad += dp->r;
 			if (div != NULL)
 				dp++;
-			else {
+			else
 				addcolor(acol, dp->v);
-				arad += dp->r;
-			}
 		}
-	if (ns > 0) {			/* perform super-sampling */
+	if (ns > 0 && arad > FTINY && ndivs/arad < minarad)
+		ns = 0;			/* close enough */
+	else if (ns > 0) {		/* else perform super-sampling */
 		comperrs(div, &hemi);			/* compute errors */
 		qsort(div, ndivs, sizeof(AMBSAMP), ambcmp);	/* sort divs */
 						/* super-sample */
@@ -170,6 +169,7 @@ FVECT  pg, dg;
 	}
 					/* compute returned values */
 	if (div != NULL) {
+		arad = 0.0;
 		for (i = ndivs, dp = div; i-- > 0; dp++) {
 			arad += dp->r;
 			if (dp->n > 1) {
@@ -206,7 +206,7 @@ FVECT  pg, dg;
 	b = 1.0/ndivs;
 	scalecolor(acol, b);
 	if (arad <= FTINY)
-		arad = FHUGE;
+		arad = maxarad;
 	else
 		arad = (ndivs+ns)/arad;
 	if (pg != NULL) {		/* reduce radius if gradient large */
@@ -214,11 +214,17 @@ FVECT  pg, dg;
 		if (d*arad*arad > 1.0)
 			arad = 1.0/sqrt(d);
 	}
-	if (arad > maxarad)
-		arad = maxarad;
-	else if (arad < minarad)
+	if (arad < minarad) {
 		arad = minarad;
-	return(arad/sqrt(wt));
+		if (pg != NULL && d*arad*arad > 1.0) {	/* cap gradient */
+			d = 1.0/arad/sqrt(d);
+			for (i = 0; i < 3; i++)
+				pg[i] *= d;
+		}
+	}
+	if ((arad /= sqrt(wt)) > maxarad)
+		arad = maxarad;
+	return(arad);
 oopsy:
 	if (div != NULL)
 		free((char *)div);
@@ -307,10 +313,10 @@ register AMBHEMI  *hp;
 posgradient(gv, da, hp)				/* compute position gradient */
 FVECT  gv;
 AMBSAMP  *da;			/* assumes standard ordering */
-AMBHEMI  *hp;
+register AMBHEMI  *hp;
 {
 	register int  i, j;
-	double  b, d;
+	double  nextsine, lastsine, b, d;
 	double  mag0, mag1;
 	double  phi, cosp, sinp, xd, yd;
 	register AMBSAMP  *dp;
@@ -319,6 +325,7 @@ AMBHEMI  *hp;
 	for (j = 0; j < hp->np; j++) {
 		dp = da + j;
 		mag0 = mag1 = 0.0;
+		lastsine = 0.0;
 		for (i = 0; i < hp->nt; i++) {
 #ifdef  DEBUG
 			if (dp->t != i || dp->p != j)
@@ -329,24 +336,26 @@ AMBHEMI  *hp;
 			if (i > 0) {
 				d = dp[-hp->np].r;
 				if (dp[0].r > d) d = dp[0].r;
-				d *= 1.0 - (double)i/hp->nt;	/* cos(t)^2 */
+							/* sin(t)*cos(t)^2 */
+				d *= lastsine * (1.0 - (double)i/hp->nt);
 				mag0 += d*(b - bright(dp[-hp->np].v));
 			}
+			nextsine = sqrt((double)(i+1)/hp->nt);
 			if (j > 0) {
 				d = dp[-1].r;
 				if (dp[0].r > d) d = dp[0].r;
-				mag1 += d*(b - bright(dp[-1].v));
+				mag1 += d * (nextsine - lastsine) *
+						(b - bright(dp[-1].v));
 			} else {
 				d = dp[hp->np-1].r;
 				if (dp[0].r > d) d = dp[0].r;
-				mag1 += d*(b - bright(dp[hp->np-1].v));
+				mag1 += d * (nextsine - lastsine) *
+						(b - bright(dp[hp->np-1].v));
 			}
 			dp += hp->np;
+			lastsine = nextsine;
 		}
-		if (hp->nt > 1) {
-			mag0 /= (double)hp->np;
-			mag1 /= (double)hp->nt;
-		}
+		mag0 *= 2.0*PI / hp->np;
 		phi = 2.0*PI * (double)j/hp->np;
 		cosp = cos(phi); sinp = sin(phi);
 		xd += mag0*cosp - mag1*sinp;
@@ -360,7 +369,7 @@ AMBHEMI  *hp;
 dirgradient(gv, da, hp)				/* compute direction gradient */
 FVECT  gv;
 AMBSAMP  *da;			/* assumes standard ordering */
-AMBHEMI  *hp;
+register AMBHEMI  *hp;
 {
 	register int  i, j;
 	double  mag;
@@ -377,7 +386,8 @@ AMBHEMI  *hp;
 				error(CONSISTENCY,
 					"division order in dirgradient");
 #endif
-			mag += sqrt((i+.5)/hp->nt)*bright(dp->v);  /* sin(t) */
+							/* tan(t) */
+			mag += bright(dp->v)/sqrt(hp->nt/(i+.5) - 1.0);
 			dp += hp->np;
 		}
 		phi = 2.0*PI * (j+.5)/hp->np + PI/2.0;
@@ -385,5 +395,5 @@ AMBHEMI  *hp;
 		yd += mag * sin(phi);
 	}
 	for (i = 0; i < 3; i++)
-		gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])*PI/(hp->nt*hp->np);
+		gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])/(hp->nt*hp->np);
 }