--- ray/src/gen/gendaymtx.c	2020/01/07 01:42:30	2.34
+++ ray/src/gen/gendaymtx.c	2022/03/11 18:04:39	2.39
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: gendaymtx.c,v 2.34 2020/01/07 01:42:30 greg Exp $";
+static const char RCSid[] = "$Id: gendaymtx.c,v 2.39 2022/03/11 18:04:39 greg Exp $";
 #endif
 /*
  *  gendaymtx.c
@@ -89,8 +89,7 @@ static const char RCSid[] = "$Id: gendaymtx.c,v 2.34 2
 #include "color.h"
 #include "sun.h"
 
-char *progname;								/* Program name */
-char errmsg[128];							/* Error message buffer */
+char *progname;					/* Program name */
 const double DC_SolarConstantE = 1367.0;	/* Solar constant W/m^2 */
 const double DC_SolarConstantL = 127.5;		/* Solar constant klux */
 
@@ -135,7 +134,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- */
@@ -252,6 +250,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? */
 
@@ -270,17 +270,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);
 
 
@@ -306,8 +309,11 @@ 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 */
@@ -368,6 +374,19 @@ main(int argc, char *argv[])
 			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;
@@ -376,19 +395,19 @@ main(int argc, char *argv[])
 			skycolor[0] = skycolor[1] = skycolor[2] = 0;
 			grefl[0] = grefl[1] = grefl[2] = 0;
 			break;
-		case 'D':			/* output suns to file */
-			sunsfp = fopen(argv[++i], "w");
-			if (!sunsfp) {
+		case 'M':			/* send sun modifiers to file */
+			if ((modsfp = fopen(argv[++i], "w")) == NULL) {
 				fprintf(stderr, "%s: cannot open '%s' for output\n",
-						argv[0], argv[i]);
+						progname, argv[i]);
 				exit(1);
 			}
-			fixed_sun_sa = PI/360.*0.533;
-			fixed_sun_sa *= PI*fixed_sun_sa;
 			break;
 		case 's':			/* sky only (no direct) */
 			suncolor[0] = suncolor[1] = suncolor[2] = 0;
 			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;
@@ -399,7 +418,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;
@@ -417,6 +436,9 @@ main(int argc, char *argv[])
 				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",
@@ -459,8 +481,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 %s\n",
+					progname, sunsfp==stdout ? "stdout" : "file");
 		if (rotation != 0)
 			fprintf(stderr, "%s: rotating output %.0f degrees\n",
 					progname, rotation);
@@ -474,6 +500,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;
@@ -483,19 +523,23 @@ 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;
 		}
-					/* 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);
+		if (!sun_in_sky && dir > (input==1 ? 20. : 20.*WHTEFFICACY))
+			fprintf(stderr,
+				"%s: warning - unusually bright at %.1f on %d-%d\n",
+					progname, hr, mo, da);
 					/* convert measured values */
-		if (dir_is_horiz && altitude > 0.)
+		if (dir_is_horiz && altitude > FTINY)
 			dir /= sin(altitude);
 		if (input == 1) {
 			dir_irrad = dir;
@@ -506,17 +550,11 @@ main(int argc, char *argv[])
 		}
 					/* compute sky patch values */
 		ComputeSky(mtx_data+mtx_offset);
-					/* output sun if indicated */
-		if (sunsfp && (altitude > 0) & (dir_illum > 1e-4)) {
-			double	srad = dir_illum/(WHTEFFICACY * fixed_sun_sa);
-			FVECT	sv;
-			vector(sv, altitude, azimuth);
-			fprintf(sunsfp, "\nvoid light solar%d\n0\n0\n", ntsteps);
-			fprintf(sunsfp, "3 %.3e %.3e %.3e\n", srad*suncolor[0],
-					srad*suncolor[1], srad*suncolor[2]);
-			fprintf(sunsfp, "\nsolar%d source sun%d\n0\n0\n", ntsteps, ntsteps);
-			fprintf(sunsfp, "4 %.6f %.6f %.6f 0.533\n", sv[0], sv[1], sv[2]);
-		}
+					/* output sun if requested */
+		if (sunsfp && sun_in_sky)
+			OutputSun(solar_minute(julian_date,hr), 1,
+						sunsfp, modsfp);
+
 		if (avgSky < 0)		/* no matrix? */
 			continue;
 
@@ -528,6 +566,7 @@ main(int argc, char *argv[])
 		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);
@@ -614,7 +653,7 @@ alldone:
 		fprintf(stderr, "%s: done.\n", progname);
 	exit(0);
 userr:
-	fprintf(stderr, "Usage: %s [-v][-h][-A][-d|-s|-n][-D file][-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:
@@ -776,6 +815,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) */