--- ray/src/gen/gendaymtx.c	2013/01/26 00:59:08	2.7
+++ ray/src/gen/gendaymtx.c	2014/06/17 21:01:21	2.16
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: gendaymtx.c,v 2.7 2013/01/26 00:59:08 greg Exp $";
+static const char RCSid[] = "$Id: gendaymtx.c,v 2.16 2014/06/17 21:01:21 greg Exp $";
 #endif
 /*
  *  gendaymtx.c
@@ -86,6 +86,7 @@ static const char RCSid[] = "$Id: gendaymtx.c,v 2.7 20
 #include <string.h>
 #include <ctype.h>
 #include "rtmath.h"
+#include "platform.h"
 #include "color.h"
 
 char *progname;								/* Program name */
@@ -108,6 +109,7 @@ double sky_clearness;			/* Sky clearness */
 double solar_rad;			/* Solar radiance */
 double sun_zenith;			/* Sun zenith angle (radians) */
 int	input = 0;				/* Input type */
+int	output = 0;				/* Output type */
 
 extern double dmax( double, double );
 extern double CalcAirMass();
@@ -208,7 +210,7 @@ static const CategoryBounds SkyClearCat[8] =
 	{ 1.950, 2.800 },
 	{ 2.800, 4.500 },
 	{ 4.500, 6.200 },
-	{ 6.200, 12.00 }	/* Clear */
+	{ 6.200, 12.01 }	/* Clear */
 };
 
 /* Luminous efficacy model coefficients */
@@ -247,7 +249,7 @@ static const ModelCoeff DirectLumEff[8] =
 };
 
 #ifndef NSUNPATCH
-#define	NSUNPATCH	4		/* # patches to spread sun into */
+#define	NSUNPATCH	4		/* max. # patches to spread sun into */
 #endif
 
 extern int jdate(int month, int day);
@@ -260,6 +262,9 @@ extern double  s_latitude;
 extern double  s_longitude;
 extern double  s_meridian;
 
+int		nsuns = NSUNPATCH;	/* number of sun patches to use */
+double		fixed_sun_sa = -1;	/* fixed solid angle per sun? */
+
 int		verbose = 0;		/* progress reports to stderr? */
 
 int		outfmt = 'a';		/* output format */
@@ -288,14 +293,33 @@ extern int	rh_init(void);
 extern float *	resize_dmatrix(float *mtx_data, int nsteps, int npatch);
 extern void	AddDirect(float *parr);
 
+
+static const char *
+getfmtname(int fmt)
+{
+	switch (fmt) {
+	case 'a':
+		return("ascii");
+	case 'f':
+		return("float");
+	case 'd':
+		return("double");
+	}
+	return("unknown");
+}
+
+
 int
 main(int argc, char *argv[])
 {
 	char	buf[256];
+	int	doheader = 1;		/* output header? */
+	double	rotation = 0;		/* site rotation (degrees) */
 	double	elevation;		/* site elevation (meters) */
 	int	dir_is_horiz;		/* direct is meas. on horizontal? */
 	float	*mtx_data = NULL;	/* our matrix data */
 	int	ntsteps = 0;		/* number of rows in matrix */
+	int	step_alloc = 0;
 	int	last_monthly = 0;	/* month of last report */
 	int	mo, da;			/* month (1-12) and day (1-31) */
 	double	hr;			/* hour (local standard time) */
@@ -315,6 +339,9 @@ main(int argc, char *argv[])
 		case 'v':			/* verbose progress reports */
 			verbose++;
 			break;
+		case 'h':			/* turn off header */
+			doheader = 0;
+			break;
 		case 'o':			/* output format */
 			switch (argv[i][2]) {
 			case 'f':
@@ -326,6 +353,20 @@ main(int argc, char *argv[])
 				goto userr;
 			}
 			break;
+		case 'O':			/* output type */
+			switch (argv[i][2]) {
+			case '0':
+				output = 0;
+				break;
+			case '1':
+				output = 1;
+				break;
+			default:
+				goto userr;
+			}
+			if (argv[i][3])
+				goto userr;
+			break;
 		case 'm':			/* Reinhart subdivisions */
 			rhsubdiv = atoi(argv[++i]);
 			break;
@@ -344,6 +385,15 @@ main(int argc, char *argv[])
 			if (skycolor[1] <= 1e-4)
 				skycolor[0] = skycolor[1] = skycolor[2] = 1;
 			break;
+		case 'r':			/* rotate distribution */
+			if (argv[i][2] && argv[i][2] != 'z')
+				goto userr;
+			rotation = atof(argv[++i]);
+			break;
+		case '5':			/* 5-phase calculation */
+			nsuns = 1;
+			fixed_sun_sa = 6.797e-05;
+			break;
 		default:
 			goto userr;
 		}
@@ -398,6 +448,9 @@ main(int argc, char *argv[])
 				progname, s_latitude, s_longitude);
 		fprintf(stderr, "%s: %d sky patches per time step\n",
 				progname, nskypatch);
+		if (rotation != 0)
+			fprintf(stderr, "%s: rotating output %.0f degrees\n",
+					progname, rotation);
 	}
 					/* convert quantities to radians */
 	s_latitude = DegToRad(s_latitude);
@@ -408,7 +461,10 @@ main(int argc, char *argv[])
 		double		sda, sta;
 					/* make space for next time step */
 		mtx_offset = 3*nskypatch*ntsteps++;
-		mtx_data = resize_dmatrix(mtx_data, ntsteps, nskypatch);
+		if (ntsteps > step_alloc) {
+			step_alloc += (step_alloc>>1) + ntsteps + 7;
+			mtx_data = resize_dmatrix(mtx_data, step_alloc, nskypatch);
+		}
 		if (dif <= 1e-4) {
 			memset(mtx_data+mtx_offset, 0, sizeof(float)*3*nskypatch);
 			continue;
@@ -421,7 +477,7 @@ main(int argc, char *argv[])
 		sda = sdec(julian_date);
 		sta = stadj(julian_date);
 		altitude = salt(sda, hr+sta);
-		azimuth = sazi(sda, hr+sta) + PI;
+		azimuth = sazi(sda, hr+sta) + PI - DegToRad(rotation);
 					/* convert measured values */
 		if (dir_is_horiz && altitude > 0.)
 			dir /= sin(altitude);
@@ -447,12 +503,25 @@ main(int argc, char *argv[])
 			break;
 		}
 					/* write out matrix */
+	if (outfmt != 'a')
+		SET_FILE_BINARY(stdout);
 #ifdef getc_unlocked
 	flockfile(stdout);
 #endif
 	if (verbose)
 		fprintf(stderr, "%s: writing %smatrix with %d time steps...\n",
 				progname, outfmt=='a' ? "" : "binary ", ntsteps);
+	if (doheader) {
+		newheader("RADIANCE", stdout);
+		printargs(argc, argv, stdout);
+		printf("LATLONG= %.8f %.8f\n", RadToDeg(s_latitude),
+					-RadToDeg(s_longitude));
+		printf("NROWS=%d\n", nskypatch);
+		printf("NCOLS=%d\n", ntsteps);
+		printf("NCOMP=3\n");
+		fputformat(getfmtname(outfmt), stdout);
+		putchar('\n');
+	}
 					/* patches are rows (outer sort) */
 	for (i = 0; i < nskypatch; i++) {
 		mtx_offset = 3*i;
@@ -494,7 +563,7 @@ main(int argc, char *argv[])
 		fprintf(stderr, "%s: done.\n", progname);
 	exit(0);
 userr:
-	fprintf(stderr, "Usage: %s [-v][-d|-s][-m N][-g r g b][-c r g b][-o{f|d}] [tape.wea]\n",
+	fprintf(stderr, "Usage: %s [-v][-h][-d|-s][-r deg][-m N][-g r g b][-c r g b][-o{f|d}][-O{0|1}] [tape.wea]\n",
 			progname);
 	exit(1);
 fmterr:
@@ -547,6 +616,14 @@ ComputeSky(float *parr)
 		sky_brightness = CalcSkyBrightness();
 		sky_clearness =  CalcSkyClearness();
 
+		/* Limit sky clearness */
+		if (sky_clearness > 11.9)
+			sky_clearness = 11.9;
+
+		/* Limit sky brightness */
+		if (sky_brightness < 0.01)
+			sky_brightness = 0.01;
+
 		/* Calculate illuminance */
 		index = GetCategoryIndex();
 		diff_illum = diff_irrad * CalcDiffuseIllumRatio(index);
@@ -558,6 +635,11 @@ ComputeSky(float *parr)
 		index = CalcSkyParamFromIllum();
 	}
 
+	if (output == 1) {			/* hack for solar radiance */
+		diff_illum = diff_irrad * WHTEFFICACY;
+		dir_illum = dir_irrad * WHTEFFICACY;
+	}
+
 	if (bright(skycolor) <= 1e-4) {			/* 0 sky component? */
 		memset(parr, 0, sizeof(float)*3*nskypatch);
 		return;
@@ -578,6 +660,12 @@ ComputeSky(float *parr)
 	/* Calculate relative horizontal illuminance */
 	norm_diff_illum = CalcRelHorzIllum(parr);
 
+	/* Check for zero sky -- make uniform in that case */
+	if (norm_diff_illum <= FTINY) {
+		for (i = 1; i < nskypatch; i++)
+			setcolor(parr+3*i, 1., 1., 1.);
+		norm_diff_illum = PI;
+	}
 	/* Normalization coefficient */
 	norm_diff_illum = diff_illum / norm_diff_illum;
 
@@ -600,8 +688,12 @@ AddDirect(float *parr)
 
 	if (dir_illum <= 1e-4 || bright(suncolor) <= 1e-4)
 		return;
-					/* identify NSUNPATCH closest patches */
-	for (i = NSUNPATCH; i--; )
+					/* identify nsuns closest patches */
+	if (nsuns > NSUNPATCH)
+		nsuns = NSUNPATCH;
+	else if (nsuns <= 0)
+		nsuns = 1;
+	for (i = nsuns; i--; )
 		near_dprod[i] = -1.;
 	vector(svec, altitude, azimuth);
 	for (p = 1; p < nskypatch; p++) {
@@ -609,9 +701,9 @@ AddDirect(float *parr)
 		double	dprod;
 		rh_vector(pvec, p);
 		dprod = DOT(pvec, svec);
-		for (i = 0; i < NSUNPATCH; i++)
+		for (i = 0; i < nsuns; i++)
 			if (dprod > near_dprod[i]) {
-				for (j = NSUNPATCH; --j > i; ) {
+				for (j = nsuns; --j > i; ) {
 					near_dprod[j] = near_dprod[j-1];
 					near_patch[j] = near_patch[j-1];
 				}
@@ -621,13 +713,15 @@ AddDirect(float *parr)
 			}
 	}
 	wtot = 0;			/* weight by proximity */
-	for (i = NSUNPATCH; i--; )
+	for (i = nsuns; i--; )
 		wtot += wta[i] = 1./(1.002 - near_dprod[i]);
 					/* add to nearest patch radiances */
-	for (i = NSUNPATCH; i--; ) {
+	for (i = nsuns; i--; ) {
 		float	*pdest = parr + 3*near_patch[i];
-		float	val_add = wta[i] * dir_illum /
-				(WHTEFFICACY * wtot * rh_dom[near_patch[i]]);
+		float	val_add = wta[i] * dir_illum / (WHTEFFICACY * wtot);
+
+		val_add /= (fixed_sun_sa > 0)	? fixed_sun_sa 
+						: rh_dom[near_patch[i]] ;
 		*pdest++ += val_add*suncolor[0];
 		*pdest++ += val_add*suncolor[1];
 		*pdest++ += val_add*suncolor[2];
@@ -771,7 +865,7 @@ double CalcSkyClearness()
 	double sz_cubed;	/* Sun zenith angle cubed */
 
 	/* Calculate sun zenith angle cubed */
-	sz_cubed = pow(sun_zenith, 3.0);
+	sz_cubed = sun_zenith*sun_zenith*sun_zenith;
 
 	return ((diff_irrad + dir_irrad) / diff_irrad + 1.041 *
 			sz_cubed) / (1.0 + 1.041 * sz_cubed);
@@ -802,7 +896,7 @@ double CalcDiffuseIrradiance()
 double CalcDirectIrradiance()
 {
 	return CalcDiffuseIrradiance() * ((sky_clearness - 1.0) * (1 + 1.041
-			* pow(sun_zenith, 3.0)));
+			* sun_zenith*sun_zenith*sun_zenith));
 }
 
 /* Calculate sky brightness and clearness from illuminance values */
@@ -828,7 +922,7 @@ int CalcSkyParamFromIllum()
 		sky_clearness = 12.0;
 
 	/* Limit sky brightness */
-	if (sky_brightness < 0.05)
+	if (sky_brightness < 0.01)
 			sky_brightness = 0.01; 
 
 	while (((fabs(diff_irrad - test1) > 10.0) ||
@@ -852,7 +946,7 @@ int CalcSkyParamFromIllum()
 			sky_clearness = 12.0;
 	
 		/* Limit sky brightness */
-		if (sky_brightness < 0.05)
+		if (sky_brightness < 0.01)
 			sky_brightness = 0.01; 
 	}