--- ray/src/gen/gendaymtx.c	2013/01/18 19:56:03	2.2
+++ ray/src/gen/gendaymtx.c	2013/01/20 01:16:15	2.4
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: gendaymtx.c,v 2.2 2013/01/18 19:56:03 greg Exp $";
+static const char RCSid[] = "$Id: gendaymtx.c,v 2.4 2013/01/20 01:16:15 greg Exp $";
 #endif
 /*
  *  gendaymtx.c
@@ -246,6 +246,10 @@ static const ModelCoeff DirectLumEff[8] =
 	{ 101.18,  1.58, -1.10,  -8.29 }
 };
 
+#ifndef NSUNPATCH
+#define	NSUNPATCH	4		/* # patches to spread sun into */
+#endif
+
 extern int jdate(int month, int day);
 extern double stadj(int  jd);
 extern double sdec(int  jd);
@@ -256,18 +260,16 @@ extern double  s_latitude;
 extern double  s_longitude;
 extern double  s_meridian;
 
-double		grefl = 0.2;		/* diffuse ground reflectance */
-
 int		verbose = 0;		/* progress reports to stderr? */
 
 int		outfmt = 'a';		/* output format */
 
 int		rhsubdiv = 1;		/* Reinhart sky subdivisions */
 
-float		skycolor[3] = {.96, 1.004, 1.118};	/* sky coloration */
+COLOR		skycolor = {.96, 1.004, 1.118};	/* sky coloration */
+COLOR		suncolor = {1., 1., 1.};	/* sun color */
+COLOR		grefl = {.2, .2, .2};		/* ground reflectance */
 
-int		do_sun = 1;		/* output direct solar contribution? */
-
 int		nskypatch;		/* number of Reinhart patches */
 float		*rh_palt;		/* sky patch altitudes (radians) */
 float		*rh_pazi;		/* sky patch azimuths (radians) */
@@ -305,13 +307,15 @@ main(int argc, char *argv[])
 					/* get options */
 	for (i = 1; i < argc && argv[i][0] == '-'; i++)
 		switch (argv[i][1]) {
-		case 'g':
-			grefl = atof(argv[++i]);
+		case 'g':			/* ground reflectance */
+			grefl[0] = atof(argv[++i]);
+			grefl[1] = atof(argv[++i]);
+			grefl[2] = atof(argv[++i]);
 			break;
-		case 'v':
+		case 'v':			/* verbose progress reports */
 			verbose++;
 			break;
-		case 'o':
+		case 'o':			/* output format */
 			switch (argv[i][2]) {
 			case 'f':
 			case 'd':
@@ -322,20 +326,21 @@ main(int argc, char *argv[])
 				goto userr;
 			}
 			break;
-		case 'm':
+		case 'm':			/* Reinhart subdivisions */
 			rhsubdiv = atoi(argv[++i]);
 			break;
-		case 'c':
+		case 'c':			/* sky color */
 			skycolor[0] = atof(argv[++i]);
 			skycolor[1] = atof(argv[++i]);
 			skycolor[2] = atof(argv[++i]);
 			break;
-		case 'd':
-			do_sun = 1;
+		case 'd':			/* solar (direct) only */
 			skycolor[0] = skycolor[1] = skycolor[2] = 0;
+			if (suncolor[1] <= 1e-4)
+				suncolor[0] = suncolor[1] = suncolor[2] = 1;
 			break;
-		case 's':
-			do_sun = 0;
+		case 's':			/* sky only (no direct) */
+			suncolor[0] = suncolor[1] = suncolor[2] = 0;
 			if (skycolor[1] <= 1e-4)
 				skycolor[0] = skycolor[1] = skycolor[2] = 1;
 			break;
@@ -429,8 +434,7 @@ main(int argc, char *argv[])
 		}
 					/* compute sky patch values */
 		ComputeSky(mtx_data+mtx_offset);
-		if (do_sun)
-			AddDirect(mtx_data+mtx_offset);
+		AddDirect(mtx_data+mtx_offset);
 	}
 					/* check for junk at end */
 	while ((i = fgetc(stdin)) != EOF)
@@ -455,7 +459,7 @@ main(int argc, char *argv[])
 		switch (outfmt) {
 		case 'a':
 			for (j = 0; j < ntsteps; j++) {
-				printf("%.3e %.3e %.3e\n", mtx_data[mtx_offset],
+				printf("%.3g %.3g %.3g\n", mtx_data[mtx_offset],
 						mtx_data[mtx_offset+1],
 						mtx_data[mtx_offset+2]);
 				mtx_offset += 3*nskypatch;
@@ -558,10 +562,11 @@ ComputeSky(float *parr)
 		return;
 	}
 	/* Compute ground radiance (include solar contribution if any) */
-	parr[0] = diff_illum * (1./PI/WHTEFFICACY);
+	parr[0] = diff_illum;
 	if (altitude > 0)
-		parr[0] += dir_illum * sin(altitude) * (1./PI/WHTEFFICACY);
-	parr[2] = parr[1] = parr[0];
+		parr[0] += dir_illum * sin(altitude);
+	parr[2] = parr[1] = parr[0] *= (1./PI/WHTEFFICACY);
+	multcolor(parr, grefl);
 
 	/* Calculate Perez sky model parameters */
 	CalcPerezParam(sun_zenith, sky_clearness, sky_brightness, index);
@@ -593,15 +598,15 @@ void
 AddDirect(float *parr)
 {
 	FVECT	svec;
-	double	near_dprod[4];
-	int	near_patch[4];
-	double	wta[4], wtot;
+	double	near_dprod[NSUNPATCH];
+	int	near_patch[NSUNPATCH];
+	double	wta[NSUNPATCH], wtot;
 	int	i, j, p;
 
-	if (!do_sun || dir_illum < 1e-4)
+	if (dir_illum <= 1e-4 || bright(suncolor) <= 1e-4)
 		return;
-					/* identify 4 closest patches */
-	for (i = 4; i--; )
+					/* identify NSUNPATCH closest patches */
+	for (i = NSUNPATCH; i--; )
 		near_dprod[i] = -1.;
 	vector(svec, altitude, azimuth);
 	for (p = 1; p < nskypatch; p++) {
@@ -609,9 +614,9 @@ AddDirect(float *parr)
 		double	dprod;
 		rh_vector(pvec, p);
 		dprod = DOT(pvec, svec);
-		for (i = 0; i < 4; i++)
+		for (i = 0; i < NSUNPATCH; i++)
 			if (dprod > near_dprod[i]) {
-				for (j = 4; --j > i; ) {
+				for (j = NSUNPATCH; --j > i; ) {
 					near_dprod[j] = near_dprod[j-1];
 					near_patch[j] = near_patch[j-1];
 				}
@@ -621,16 +626,16 @@ AddDirect(float *parr)
 			}
 	}
 	wtot = 0;			/* weight by proximity */
-	for (i = 4; i--; )
+	for (i = NSUNPATCH; i--; )
 		wtot += wta[i] = 1./(1.002 - near_dprod[i]);
 					/* add to nearest patch radiances */
-	for (i = 4; i--; ) {
+	for (i = NSUNPATCH; i--; ) {
 		float	*pdest = parr + 3*near_patch[i];
 		float	val_add = wta[i] * dir_illum /
 				(WHTEFFICACY * wtot * rh_dom[near_patch[i]]);
-		*pdest++ += val_add;
-		*pdest++ += val_add;
-		*pdest++ += val_add;
+		*pdest++ += val_add*suncolor[0];
+		*pdest++ += val_add*suncolor[1];
+		*pdest++ += val_add*suncolor[2];
 	}
 }
 
@@ -665,7 +670,8 @@ rh_init(void)
 	for (i = 0; i < NROW*rhsubdiv; i++) {
 		const float	ralt = alpha*(i + .5);
 		const int	ninrow = tnaz[i/rhsubdiv]*rhsubdiv;
-		const float	dom = (sin(alpha*(i+1)) - sin(alpha*i))/ninrow;
+		const float	dom = 2.*PI*(sin(alpha*(i+1)) - sin(alpha*i)) /
+						(double)ninrow;
 		for (j = 0; j < ninrow; j++) {
 			rh_palt[p] = ralt;
 			rh_pazi[p] = 2.*PI * j / (double)ninrow;