--- ray/src/gen/gendaymtx.c	2013/01/20 01:16:49	2.5
+++ ray/src/gen/gendaymtx.c	2013/02/05 06:00:19	2.8
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: gendaymtx.c,v 2.5 2013/01/20 01:16:49 greg Exp $";
+static const char RCSid[] = "$Id: gendaymtx.c,v 2.8 2013/02/05 06:00:19 greg Exp $";
 #endif
 /*
  *  gendaymtx.c
@@ -292,6 +292,7 @@ int
 main(int argc, char *argv[])
 {
 	char	buf[256];
+	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 */
@@ -344,6 +345,11 @@ 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;
 		default:
 			goto userr;
 		}
@@ -398,6 +404,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);
@@ -421,7 +430,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);
@@ -494,7 +503,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][-d|-s][-r deg][-m N][-g r g b][-c r g b][-o{f|d}] [tape.wea]\n",
 			progname);
 	exit(1);
 fmterr:
@@ -515,19 +524,20 @@ ComputeSky(float *parr)
 {
 	int index;			/* Category index */
 	double norm_diff_illum;		/* Normalized diffuse illuimnance */
-	double zlumin;			/* Zenith luminance */
 	int i;
 	
 	/* Calculate atmospheric precipitable water content */
 	apwc = CalcPrecipWater(dew_point);
 
-	/* Limit solar altitude to keep circumsolar off zenith */
-	if (altitude > DegToRad(87.0))
-		altitude = DegToRad(87.0);
+	/* Calculate sun zenith angle (don't let it dip below horizon) */
+	/* Also limit minimum angle to keep circumsolar off zenith */
+	if (altitude <= 0.0)
+		sun_zenith = DegToRad(90.0);
+	else if (altitude >= DegToRad(87.0))
+		sun_zenith = DegToRad(3.0);
+	else
+		sun_zenith = DegToRad(90.0) - altitude;
 
-	/* Calculate sun zenith angle */
-	sun_zenith = DegToRad(90.0) - altitude;
-
 	/* Compute the inputs for the calculation of the sky distribution */
 	
 	if (input == 0)					/* XXX never used */
@@ -580,15 +590,9 @@ ComputeSky(float *parr)
 	/* Normalization coefficient */
 	norm_diff_illum = diff_illum / norm_diff_illum;
 
-	/* Calculate relative zenith luminance */
-	zlumin = CalcRelLuminance(sun_zenith, 0.0);
-
-	/* Calculate absolute zenith illuminance */
-	zlumin *= norm_diff_illum;
-
 	/* Apply to sky patches to get absolute radiance values */
 	for (i = 1; i < nskypatch; i++) {
-		scalecolor(parr+3*i, zlumin*(1./WHTEFFICACY));
+		scalecolor(parr+3*i, norm_diff_illum*(1./WHTEFFICACY));
 		multcolor(parr+3*i, skycolor);
 	}
 }
@@ -943,9 +947,9 @@ double CalcRelHorzIllum( float *parr )
 	double rh_illum = 0.0;	/* Relative horizontal illuminance */
 
 	for (i = 1; i < nskypatch; i++)
-		rh_illum += parr[3*i+1] * rh_cos(i);
+		rh_illum += parr[3*i+1] * rh_cos(i) * rh_dom[i];
 
-	return rh_illum * (2.0 * PI / (nskypatch-1));
+	return rh_illum;
 }
 
 /* Calculate earth orbit eccentricity correction factor */