--- ray/src/cv/bsdfmesh.c	2014/08/21 10:33:48	2.30
+++ ray/src/cv/bsdfmesh.c	2014/08/22 05:38:44	2.33
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfmesh.c,v 2.30 2014/08/21 10:33:48 greg Exp $";
+static const char RCSid[] = "$Id: bsdfmesh.c,v 2.33 2014/08/22 05:38:44 greg Exp $";
 #endif
 /*
  * Create BSDF advection mesh from radial basis functions.
@@ -51,7 +51,7 @@ eval_DSFsurround(const RBFNODE *rbf, const FVECT outve
 		if (i) spinvector(tvec, tvec, outvec, phinc);
 		if (tvec[2] > 0 ^ output_orient > 0)
 			continue;
-		sum += eval_rbfrep(rbf, tvec) * output_orient*tvec[2];
+		sum += eval_rbfrep(rbf, tvec) * COSF(tvec[2]);
 		++n;
 	}
 	if (n < 2)				/* should never happen! */
@@ -64,26 +64,33 @@ static double
 est_DSFrad(const RBFNODE *rbf, const FVECT outvec)
 {
 	const double	rad_epsilon = 0.03;
-	const double	DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec)
-						* output_orient*outvec[2];
+	const double	DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec) *
+							COSF(outvec[2]);
 	double		inside_rad = rad_epsilon;
 	double		outside_rad = 0.5;
 	double		DSFinside = eval_DSFsurround(rbf, outvec, inside_rad);
 	double		DSFoutside = eval_DSFsurround(rbf, outvec, outside_rad);
 #define	interp_rad	inside_rad + (outside_rad-inside_rad) * \
 				(DSFtarget-DSFinside) / (DSFoutside-DSFinside)
-						/* interpolation search */
-	while (outside_rad-inside_rad > rad_epsilon) {
+						/* Newton's method (sort of) */
+	do {
 		double	test_rad = interp_rad;
-		double	DSFtest = eval_DSFsurround(rbf, outvec, test_rad);
-		if (DSFtarget < DSFtest) {
+		double	DSFtest;
+		if (test_rad >= outside_rad)
+			return(test_rad);
+		if (test_rad <= inside_rad)
+			return(test_rad*(test_rad>0));
+		DSFtest = eval_DSFsurround(rbf, outvec, test_rad);
+		if (DSFtest > DSFtarget) {
 			inside_rad = test_rad;
 			DSFinside = DSFtest;
 		} else {
 			outside_rad = test_rad;
 			DSFoutside = DSFtest;
 		}
-	}
+		if (DSFoutside >= DSFinside)
+			return(test_rad);
+	} while (outside_rad-inside_rad > rad_epsilon);
 	return(interp_rad);
 #undef interp_rad
 }
@@ -100,7 +107,7 @@ comp_bsdf_spec(void)
 
 	if (dsf_list == NULL) {
 		bsdf_spec_peak = 0;
-		bsdf_spec_crad = 0;
+		bsdf_spec_rad = 0;
 		return;
 	}
 	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
@@ -112,7 +119,7 @@ comp_bsdf_spec(void)
 		++n;
 	}
 	bsdf_spec_peak = peak_sum/(double)n;
-	bsdf_spec_crad = ANG2R( rad_sum/(double)n );
+	bsdf_spec_rad = rad_sum/(double)n;
 }
 
 /* Create a new migration holder (sharing memory for multiprocessing) */