--- ray/src/gen/gendaymtx.c	2013/01/25 23:17:31	2.6
+++ ray/src/gen/gendaymtx.c	2013/02/09 00:20:24	2.9
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: gendaymtx.c,v 2.6 2013/01/25 23:17:31 greg Exp $";
+static const char RCSid[] = "$Id: gendaymtx.c,v 2.9 2013/02/09 00:20:24 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,7 +524,6 @@ 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 */
@@ -548,6 +556,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);
@@ -582,15 +598,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);
 	}
 }
@@ -835,7 +845,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) ||
@@ -859,7 +869,7 @@ int CalcSkyParamFromIllum()
 			sky_clearness = 12.0;
 	
 		/* Limit sky brightness */
-		if (sky_brightness < 0.05)
+		if (sky_brightness < 0.01)
 			sky_brightness = 0.01; 
 	}
 
@@ -945,9 +955,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 */