--- ray/src/gen/gendaymtx.c	2019/11/07 23:15:07	2.30
+++ ray/src/gen/gendaymtx.c	2020/08/15 03:28:56	2.37
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: gendaymtx.c,v 2.30 2019/11/07 23:15:07 greg Exp $";
+static const char RCSid[] = "$Id: gendaymtx.c,v 2.37 2020/08/15 03:28:56 greg Exp $";
 #endif
 /*
  *  gendaymtx.c
@@ -81,9 +81,7 @@ static const char RCSid[] = "$Id: gendaymtx.c,v 2.30 2
 
 /* Include files */
 #define	_USE_MATH_DEFINES
-#include <stdio.h>
 #include <stdlib.h>
-#include <string.h>
 #include <ctype.h>
 #include "platform.h"
 #include "rtmath.h"
@@ -137,7 +135,6 @@ extern void ComputeSky( float *parr );
 /* Radiuans into degrees */
 #define RadToDeg(rad)	((rad)*(180./PI))
 
-
 /* Perez sky model coefficients */
 
 /* Reference:	Perez, R., R. Seals, and J. Michalsky, 1993. "All- */
@@ -254,6 +251,8 @@ static const ModelCoeff DirectLumEff[8] =
 #define	NSUNPATCH	4		/* max. # patches to spread sun into */
 #endif
 
+#define	SUN_ANG_DEG	0.533		/* sun full-angle in degrees */
+
 int		nsuns = NSUNPATCH;	/* number of sun patches to use */
 double		fixed_sun_sa = -1;	/* fixed solid angle per sun? */
 
@@ -272,17 +271,20 @@ float		*rh_palt;		/* sky patch altitudes (radians) */
 float		*rh_pazi;		/* sky patch azimuths (radians) */
 float		*rh_dom;		/* sky patch solid angle (sr) */
 
-#define		vector(v,alt,azi)	(	(v)[1] = tcos(alt), \
-						(v)[0] = (v)[1]*tsin(azi), \
-						(v)[1] *= tcos(azi), \
-						(v)[2] = tsin(alt) )
+#define		vector(v,alt,azi)	(	(v)[1] = cos(alt), \
+						(v)[0] = (v)[1]*sin(azi), \
+						(v)[1] *= cos(azi), \
+						(v)[2] = sin(alt) )
 
 #define		rh_vector(v,i)		vector(v,rh_palt[i],rh_pazi[i])
 
 #define		rh_cos(i)		tsin(rh_palt[i])
 
+#define		solar_minute(jd,hr)	((24*60)*((jd)-1)+(int)((hr)*60.+.5))
+
 extern int	rh_init(void);
 extern float *	resize_dmatrix(float *mtx_data, int nsteps, int npatch);
+extern void	OutputSun(int id, int goodsun, FILE *fp, FILE *mfp);
 extern void	AddDirect(float *parr);
 
 
@@ -308,14 +310,17 @@ main(int argc, char *argv[])
 	int	doheader = 1;		/* output header? */
 	double	rotation = 0;		/* site rotation (degrees) */
 	double	elevation;		/* site elevation (meters) */
+	int	leap_day = 0;		/* add leap day? */
+	int	sun_hours_only = 0;	/* only output sun hours? */
 	int	dir_is_horiz;		/* direct is meas. on horizontal? */
+	FILE	*sunsfp = NULL;		/* output file for individual suns */
+	FILE	*modsfp = NULL;		/* modifier output file */
 	float	*mtx_data = NULL;	/* our matrix data */
 	int	avgSky = 0;		/* compute average sky r.t. matrix? */
 	int	ntsteps = 0;		/* number of time steps */
 	int	tstorage = 0;		/* number of allocated time steps */
 	int	nstored = 0;		/* number of time steps in matrix */
 	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 */
@@ -366,25 +371,44 @@ 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':			/* output suns to file */
+			if (strcmp(argv[++i], "-")) {
+				sunsfp = fopen(argv[i], "w");
+				if (sunsfp == NULL) {
+					fprintf(stderr,
+					"%s: cannot open '%s' for output\n",
+							progname, argv[i]);
+					exit(1);
+				}
+				break;		/* still may output matrix */
+			}
+			sunsfp = stdout;	/* sending to stdout, so... */
+			/* fall through */
+		case 'n':			/* no matrix output */
+			avgSky = -1;
+			rhsubdiv = 1;
+			/* fall through */
 		case 'd':			/* solar (direct) only */
 			skycolor[0] = skycolor[1] = skycolor[2] = 0;
-			if (suncolor[1] <= 1e-4) {
-				inconsistent = 1;
-				suncolor[0] = suncolor[1] = suncolor[2] = 1;
+			grefl[0] = grefl[1] = grefl[2] = 0;
+			break;
+		case 'M':			/* send sun modifiers to file */
+			if ((modsfp = fopen(argv[++i], "w")) == NULL) {
+				fprintf(stderr, "%s: cannot open '%s' for output\n",
+						progname, argv[i]);
+				exit(1);
 			}
 			break;
 		case 's':			/* sky only (no direct) */
 			suncolor[0] = suncolor[1] = suncolor[2] = 0;
-			if (skycolor[1] <= 1e-4) {
-				inconsistent = 1;
-				skycolor[0] = skycolor[1] = skycolor[2] = 1;
-			}
 			break;
+		case 'u':			/* solar hours only */
+			sun_hours_only = 1;
+			break;
 		case 'r':			/* rotate distribution */
 			if (argv[i][2] && argv[i][2] != 'z')
 				goto userr;
@@ -395,7 +419,7 @@ main(int argc, char *argv[])
 			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]);
+						progname);
 				exit(1);
 			}
 			fixed_sun_sa *= fixed_sun_sa*PI;
@@ -408,14 +432,14 @@ main(int argc, char *argv[])
 		}
 	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]);
 		exit(1);
 	}
+	if ((modsfp != NULL) & (sunsfp == NULL))
+		fprintf(stderr, "%s: warning -M output will be empty without -D\n",
+				progname);
 	if (verbose) {
 		if (i == argc-1)
 			fprintf(stderr, "%s: reading weather tape '%s'\n",
@@ -458,8 +482,12 @@ main(int argc, char *argv[])
 		fprintf(stderr, "%s: location '%s'\n", progname, buf);
 		fprintf(stderr, "%s: (lat,long)=(%.1f,%.1f) degrees north, west\n",
 				progname, s_latitude, s_longitude);
-		fprintf(stderr, "%s: %d sky patches per time step\n",
-				progname, nskypatch);
+		if (avgSky >= 0)
+			fprintf(stderr, "%s: %d sky patches\n",
+					progname, nskypatch);
+		if (sunsfp)
+			fprintf(stderr, "%s: outputting suns to file\n",
+					progname);
 		if (rotation != 0)
 			fprintf(stderr, "%s: rotating output %.0f degrees\n",
 					progname, rotation);
@@ -473,6 +501,20 @@ main(int argc, char *argv[])
 					/* process each time step in tape */
 	while (scanf("%d %d %lf %lf %lf\n", &mo, &da, &hr, &dir, &dif) == 5) {
 		double		sda, sta;
+		int		sun_in_sky;
+					/* compute solar position */
+		if ((mo == 2) & (da == 29)) {
+			julian_date = 60;
+			leap_day = 1;
+		} else
+			julian_date = jdate(mo, da) + leap_day;
+		sda = sdec(julian_date);
+		sta = stadj(julian_date);
+		altitude = salt(sda, hr+sta);
+		sun_in_sky = (altitude > -DegToRad(SUN_ANG_DEG/2.));
+		if (sun_hours_only && !sun_in_sky)
+			continue;	/* skipping nighttime points */
+		azimuth = sazi(sda, hr+sta) + PI - DegToRad(rotation);
 
 		mtx_offset = 3*nskypatch*nstored;
 		nstored += !avgSky | !nstored;
@@ -482,20 +524,17 @@ main(int argc, char *argv[])
 			mtx_data = resize_dmatrix(mtx_data, tstorage, nskypatch);
 		}
 		ntsteps++;		/* keep count of time steps */
-		if (dif <= 1e-4) {
+
+		if (dir+dif <= 1e-4) {	/* effectively nighttime? */
 			if (!avgSky | !mtx_offset)
-				memset(mtx_data+mtx_offset, 0, sizeof(float)*3*nskypatch);
+				memset(mtx_data+mtx_offset, 0,
+						sizeof(float)*3*nskypatch);
+					/* output black sun? */
+			if (sunsfp && sun_in_sky)
+				OutputSun(solar_minute(julian_date,hr), 0,
+							sunsfp, modsfp);
 			continue;
 		}
-		if (verbose && mo != last_monthly)
-			fprintf(stderr, "%s: stepping through month %d...\n",
-						progname, last_monthly=mo);
-					/* compute solar position */
-		julian_date = jdate(mo, da);
-		sda = sdec(julian_date);
-		sta = stadj(julian_date);
-		altitude = salt(sda, hr+sta);
-		azimuth = sazi(sda, hr+sta) + PI - DegToRad(rotation);
 					/* convert measured values */
 		if (dir_is_horiz && altitude > 0.)
 			dir /= sin(altitude);
@@ -508,11 +547,28 @@ main(int argc, char *argv[])
 		}
 					/* compute sky patch values */
 		ComputeSky(mtx_data+mtx_offset);
+					/* output sun if requested */
+		if (sunsfp && sun_in_sky)
+			OutputSun(solar_minute(julian_date,hr), 1,
+						sunsfp, modsfp);
+
+		if (avgSky < 0)		/* no matrix? */
+			continue;
+
 		AddDirect(mtx_data+mtx_offset);
 					/* update cumulative sky? */
 		for (i = 3*nskypatch*(avgSky&(ntsteps>1)); i--; )
 			mtx_data[i] += mtx_data[mtx_offset+i];
+					/* monthly reporting */
+		if (verbose && mo != last_monthly)
+			fprintf(stderr, "%s: stepping through month %d...\n",
+						progname, last_monthly=mo);
+					/* note whether leap-day was given */
 	}
+	if (!ntsteps) {
+		fprintf(stderr, "%s: no valid time steps on input\n", progname);
+		exit(1);
+	}
 					/* check for junk at end */
 	while ((i = fgetc(stdin)) != EOF)
 		if (!isspace(i)) {
@@ -523,10 +579,10 @@ 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);
-	}
+
+	if (avgSky < 0)			/* no matrix output? */
+		goto alldone;
+
 	dif = 1./(double)ntsteps;	/* average sky? */
 	for (i = 3*nskypatch*(avgSky&(ntsteps>1)); i--; )
 		mtx_data[i] *= dif;
@@ -587,17 +643,18 @@ main(int argc, char *argv[])
 		if (ferror(stdout))
 			goto writerr;
 	}
-	if (fflush(stdout) == EOF)
+alldone:
+	if (fflush(NULL) == EOF)
 		goto writerr;
 	if (verbose)
 		fprintf(stderr, "%s: done.\n", progname);
 	exit(0);
 userr:
-	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",
+	fprintf(stderr, "Usage: %s [-v][-h][-A][-d|-s|-n][-u][-D file [-M modfile]][-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:
-	fprintf(stderr, "%s: input weather tape format error\n", progname);
+	fprintf(stderr, "%s: weather tape format error in header\n", progname);
 	exit(1);
 writerr:
 	fprintf(stderr, "%s: write error on output\n", progname);
@@ -669,11 +726,6 @@ ComputeSky(float *parr)
 		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;
-	}
 	/* Compute ground radiance (include solar contribution if any) */
 	parr[0] = diff_illum;
 	if (altitude > 0)
@@ -681,6 +733,10 @@ ComputeSky(float *parr)
 	parr[2] = parr[1] = parr[0] *= (1./PI/WHTEFFICACY);
 	multcolor(parr, grefl);
 
+	if (bright(skycolor) <= 1e-4) {			/* 0 sky component? */
+		memset(parr+3, 0, sizeof(float)*3*(nskypatch-1));
+		return;
+	}
 	/* Calculate Perez sky model parameters */
 	CalcPerezParam(sun_zenith, sky_clearness, sky_brightness, index);
 
@@ -756,6 +812,26 @@ AddDirect(float *parr)
 		*pdest++ += val_add*suncolor[1];
 		*pdest++ += val_add*suncolor[2];
 	}
+}
+
+/* Output a sun to indicated file if appropriate for this time step */
+void
+OutputSun(int id, int goodsun, FILE *fp, FILE *mfp)
+{
+	double	srad;
+	FVECT	sv;
+
+	srad = DegToRad(SUN_ANG_DEG/2.);
+	srad = goodsun ? dir_illum/(WHTEFFICACY * PI*srad*srad) : 0;
+	vector(sv, altitude, azimuth);
+	fprintf(fp, "\nvoid light solar%d\n0\n0\n", id);
+	fprintf(fp, "3 %.3e %.3e %.3e\n", srad*suncolor[0],
+			srad*suncolor[1], srad*suncolor[2]);
+	fprintf(fp, "\nsolar%d source sun%d\n0\n0\n", id, id);
+	fprintf(fp, "4 %.6f %.6f %.6f %.4f\n", sv[0], sv[1], sv[2], SUN_ANG_DEG);
+	
+	if (mfp != NULL)		/* saving modifier IDs? */
+		fprintf(mfp, "solar%d\n", id);
 }
 
 /* Initialize Reinhart sky patch positions (GW) */