--- ray/src/gen/gendaymtx.c	2015/08/01 23:27:04	2.20
+++ ray/src/gen/gendaymtx.c	2019/06/25 00:09:45	2.27
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: gendaymtx.c,v 2.20 2015/08/01 23:27:04 greg Exp $";
+static const char RCSid[] = "$Id: gendaymtx.c,v 2.27 2019/06/25 00:09:45 greg Exp $";
 #endif
 /*
  *  gendaymtx.c
@@ -86,6 +86,7 @@ static const char RCSid[] = "$Id: gendaymtx.c,v 2.20 2
 #include <string.h>
 #include <ctype.h>
 #include "rtmath.h"
+#include "rtio.h"
 #include "resolu.h"
 #include "platform.h"
 #include "color.h"
@@ -320,9 +321,12 @@ main(int argc, char *argv[])
 	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	ntsteps = 0;		/* number of time steps */
+	int	avgSky = 0;		/* compute average sky r.t. matrix? */
+	int	nrows = 0;		/* number of rows in matrix */
+	int	rows_alloc = 0;		/* number of allocated rows */
 	int	last_monthly = 0;	/* month of last report */
+	int	inconsistent = 0;	/* inconsistent options set? */
 	int	mo, da;			/* month (1-12) and day (1-31) */
 	double	hr;			/* hour (local standard time) */
 	double	dir, dif;		/* direct and diffuse values */
@@ -373,19 +377,24 @@ main(int argc, char *argv[])
 			rhsubdiv = atoi(argv[++i]);
 			break;
 		case 'c':			/* sky color */
+			inconsistent |= (skycolor[1] <= 1e-4);
 			skycolor[0] = atof(argv[++i]);
 			skycolor[1] = atof(argv[++i]);
 			skycolor[2] = atof(argv[++i]);
 			break;
 		case 'd':			/* solar (direct) only */
 			skycolor[0] = skycolor[1] = skycolor[2] = 0;
-			if (suncolor[1] <= 1e-4)
+			if (suncolor[1] <= 1e-4) {
+				inconsistent = 1;
 				suncolor[0] = suncolor[1] = suncolor[2] = 1;
+			}
 			break;
 		case 's':			/* sky only (no direct) */
 			suncolor[0] = suncolor[1] = suncolor[2] = 0;
-			if (skycolor[1] <= 1e-4)
+			if (skycolor[1] <= 1e-4) {
+				inconsistent = 1;
 				skycolor[0] = skycolor[1] = skycolor[2] = 1;
+			}
 			break;
 		case 'r':			/* rotate distribution */
 			if (argv[i][2] && argv[i][2] != 'z')
@@ -395,13 +404,24 @@ main(int argc, char *argv[])
 		case '5':			/* 5-phase calculation */
 			nsuns = 1;
 			fixed_sun_sa = PI/360.*atof(argv[++i]);
+			if (fixed_sun_sa <= 0) {
+				fprintf(stderr, "%s: missing solar disk size argument for '-5' option\n",
+						argv[0]);
+				exit(1);
+			}
 			fixed_sun_sa *= fixed_sun_sa*PI;
 			break;
+		case 'A':			/* compute average sky */
+			avgSky = 1;
+			break;
 		default:
 			goto userr;
 		}
 	if (i < argc-1)
 		goto userr;
+	if (inconsistent)
+		fprintf(stderr, "%s: WARNING: inconsistent -s, -d, -c options!\n",
+				progname);
 	if (i == argc-1 && freopen(argv[i], "r", stdin) == NULL) {
 		fprintf(stderr, "%s: cannot open '%s' for input\n",
 				progname, argv[i]);
@@ -459,17 +479,23 @@ main(int argc, char *argv[])
 	s_latitude = DegToRad(s_latitude);
 	s_longitude = DegToRad(s_longitude);
 	s_meridian = DegToRad(s_meridian);
+					/* initial allocation */
+	mtx_data = resize_dmatrix(mtx_data, rows_alloc=2, nskypatch);
 					/* process each time step in tape */
 	while (scanf("%d %d %lf %lf %lf\n", &mo, &da, &hr, &dir, &dif) == 5) {
 		double		sda, sta;
-					/* make space for next time step */
-		mtx_offset = 3*nskypatch*ntsteps++;
-		if (ntsteps > step_alloc) {
-			step_alloc += (step_alloc>>1) + ntsteps + 7;
-			mtx_data = resize_dmatrix(mtx_data, step_alloc, nskypatch);
+
+		mtx_offset = 3*nskypatch*nrows;
+		nrows += !avgSky | !nrows;
+					/* make space for next row */
+		if (nrows > rows_alloc) {
+			rows_alloc += (rows_alloc>>1) + nrows + 7;
+			mtx_data = resize_dmatrix(mtx_data, rows_alloc, nskypatch);
 		}
+		ntsteps++;		/* keep count of time steps */
 		if (dif <= 1e-4) {
-			memset(mtx_data+mtx_offset, 0, sizeof(float)*3*nskypatch);
+			if (!avgSky | !nrows)
+				memset(mtx_data+mtx_offset, 0, sizeof(float)*3*nskypatch);
 			continue;
 		}
 		if (verbose && mo != last_monthly)
@@ -494,6 +520,9 @@ main(int argc, char *argv[])
 					/* compute sky patch values */
 		ComputeSky(mtx_data+mtx_offset);
 		AddDirect(mtx_data+mtx_offset);
+					/* update cumulative sky? */
+		for (i = 3*nskypatch*(avgSky&(ntsteps>1)); i--; )
+			mtx_data[i] += mtx_data[mtx_offset+i];
 	}
 					/* check for junk at end */
 	while ((i = fgetc(stdin)) != EOF)
@@ -505,6 +534,13 @@ main(int argc, char *argv[])
 			fputs(buf, stderr); fputc('\n', stderr);
 			break;
 		}
+	if (!ntsteps) {
+		fprintf(stderr, "%s: no valid time steps on input\n", progname);
+		exit(1);
+	}
+	dif = 1./(double)ntsteps;	/* average sky? */
+	for (i = 3*nskypatch*(avgSky&(ntsteps>1)); i--; )
+		mtx_data[i] *= dif;
 					/* write out matrix */
 	if (outfmt != 'a')
 		SET_FILE_BINARY(stdout);
@@ -513,14 +549,14 @@ main(int argc, char *argv[])
 #endif
 	if (verbose)
 		fprintf(stderr, "%s: writing %smatrix with %d time steps...\n",
-				progname, outfmt=='a' ? "" : "binary ", ntsteps);
+				progname, outfmt=='a' ? "" : "binary ", nrows);
 	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("NCOLS=%d\n", nrows);
 		printf("NCOMP=3\n");
 		fputformat((char *)getfmtname(outfmt), stdout);
 		putchar('\n');
@@ -530,29 +566,29 @@ main(int argc, char *argv[])
 		mtx_offset = 3*i;
 		switch (outfmt) {
 		case 'a':
-			for (j = 0; j < ntsteps; j++) {
+			for (j = 0; j < nrows; j++) {
 				printf("%.3g %.3g %.3g\n", mtx_data[mtx_offset],
 						mtx_data[mtx_offset+1],
 						mtx_data[mtx_offset+2]);
 				mtx_offset += 3*nskypatch;
 			}
-			if (ntsteps > 1)
+			if (nrows > 1)
 				fputc('\n', stdout);
 			break;
 		case 'f':
-			for (j = 0; j < ntsteps; j++) {
-				fwrite(mtx_data+mtx_offset, sizeof(float), 3,
+			for (j = 0; j < nrows; j++) {
+				putbinary(mtx_data+mtx_offset, sizeof(float), 3,
 						stdout);
 				mtx_offset += 3*nskypatch;
 			}
 			break;
 		case 'd':
-			for (j = 0; j < ntsteps; j++) {
+			for (j = 0; j < nrows; j++) {
 				double	ment[3];
 				ment[0] = mtx_data[mtx_offset];
 				ment[1] = mtx_data[mtx_offset+1];
 				ment[2] = mtx_data[mtx_offset+2];
-				fwrite(ment, sizeof(double), 3, stdout);
+				putbinary(ment, sizeof(double), 3, stdout);
 				mtx_offset += 3*nskypatch;
 			}
 			break;
@@ -566,7 +602,7 @@ main(int argc, char *argv[])
 		fprintf(stderr, "%s: done.\n", progname);
 	exit(0);
 userr:
-	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",
+	fprintf(stderr, "Usage: %s [-v][-h][-A][-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:
@@ -938,7 +974,9 @@ int CalcSkyParamFromIllum()
 		/* Convert illuminance to irradiance */
 		index = GetCategoryIndex();
 		diff_irrad = diff_illum / CalcDiffuseIllumRatio(index);
-		dir_irrad = dir_illum / CalcDirectIllumRatio(index);
+		dir_irrad = CalcDirectIllumRatio(index);
+		if (dir_irrad > 0.1)
+			dir_irrad = dir_illum / dir_irrad;
 	
 		/* Calculate sky brightness and clearness */
 		sky_brightness = CalcSkyBrightness();