--- ray/src/cv/bsdfmesh.c	2014/08/21 10:33:48	2.30
+++ ray/src/cv/bsdfmesh.c	2017/05/16 20:41:03	2.38
@@ -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.38 2017/05/16 20:41:03 greg Exp $";
 #endif
 /*
  * Create BSDF advection mesh from radial basis functions.
@@ -7,7 +7,7 @@ static const char RCSid[] = "$Id: bsdfmesh.c,v 2.30 20
  *	G. Ward
  */
 
-#ifndef _WIN32
+#if !defined(_WIN32) && !defined(_WIN64)
 #include <unistd.h>
 #include <sys/wait.h>
 #include <sys/mman.h>
@@ -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,55 +64,95 @@ 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
 }
 
-/* Compute average BSDF peak from current DSF's */
+static int
+dbl_cmp(const void *p1, const void *p2)
+{
+	double	d1 = *(const double *)p1;
+	double	d2 = *(const double *)p2;
+
+	if (d1 > d2) return(1);
+	if (d1 < d2) return(-1);
+	return(0);
+}
+
+/* Conservative estimate of average BSDF value from current DSF's */
 static void
 comp_bsdf_spec(void)
 {
-	double		peak_sum = 0;
+	double		vmod_sum = 0;
 	double		rad_sum = 0;
 	int		n = 0;
+	double		*cost_list = NULL;
+	double		max_cost = 1.;
 	RBFNODE		*rbf;
 	FVECT		sdv;
-
-	if (dsf_list == NULL) {
-		bsdf_spec_peak = 0;
-		bsdf_spec_crad = 0;
+						/* sort by incident altitude */
+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
+		n++;
+	if (n >= 10)
+		cost_list = (double *)malloc(sizeof(double)*n);
+	if (cost_list == NULL) {
+		bsdf_spec_val = 0;
+		bsdf_spec_rad = 0;
 		return;
 	}
+	n = 0;
+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
+		cost_list[n++] = rbf->invec[2]*input_orient;
+	qsort(cost_list, n, sizeof(double), dbl_cmp);
+	max_cost = cost_list[(n+3)/4];		/* accept 25% nearest grazing */
+	free(cost_list);
+	n = 0;
 	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
+		double	this_rad, cosfact, vest;
+		if (rbf->invec[2]*input_orient > max_cost)
+			continue;
 		sdv[0] = -rbf->invec[0];
 		sdv[1] = -rbf->invec[1];
 		sdv[2] = rbf->invec[2]*(2*(input_orient==output_orient) - 1);
-		peak_sum += eval_rbfrep(rbf, sdv);
-		rad_sum += est_DSFrad(rbf, sdv);
+		cosfact = COSF(sdv[2]);
+		this_rad = est_DSFrad(rbf, sdv);
+		vest = eval_rbfrep(rbf, sdv) * cosfact *
+				(2.*M_PI) * this_rad*this_rad;
+		if (vest > rbf->vtotal)		/* don't over-estimate energy */
+			vest = rbf->vtotal;
+		vmod_sum += vest / cosfact;	/* remove cosine factor */
+		rad_sum += this_rad;
 		++n;
 	}
-	bsdf_spec_peak = peak_sum/(double)n;
-	bsdf_spec_crad = ANG2R( rad_sum/(double)n );
+	bsdf_spec_rad = rad_sum/(double)n;
+	bsdf_spec_val = vmod_sum/(2.*M_PI*n*bsdf_spec_rad*bsdf_spec_rad);
 }
 
 /* Create a new migration holder (sharing memory for multiprocessing) */
@@ -122,7 +162,7 @@ new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
 	size_t		memlen = sizeof(MIGRATION) +
 				sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
 	MIGRATION	*newmig;
-#ifdef _WIN32
+#if defined(_WIN32) || defined(_WIN64)
 	if (nprocs > 1)
 		fprintf(stderr, "%s: warning - multiprocessing not supported\n",
 				progname);
@@ -153,7 +193,7 @@ new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
 	return(mig_list = newmig);
 }
 
-#ifdef _WIN32
+#if defined(_WIN32) || defined(_WIN64)
 #define await_children(n)	(void)(n)
 #define run_subprocess()	0
 #define end_subprocess()	(void)0
@@ -451,7 +491,7 @@ check_normal_incidence(void)
 		default:
 			return;			/* else we can interpolate */
 		}
-		for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
+		for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
 			const double	d = input_orient*rbf->invec[2];
 			if (d >= 1.-2.*FTINY)
 				return;		/* seems we have normal */