--- ray/src/gen/skybright.cal	1991/11/12 17:05:15	2.1
+++ ray/src/gen/skybright.cal	2003/02/22 02:07:24	2.7
@@ -1,29 +1,38 @@
-{ SCCSid "$SunId$ LBL" }
-
+{ RCSid: $Id: skybright.cal,v 2.7 2003/02/22 02:07:24 greg Exp $ }
 {
      Sky brightness function for sunny and cloudy skies.
 
-	6/6/86
-
      Additional arguments required for calculation of skybright:
 
-	A1			- boolean, positive for cloudy sky
+	A1			- 1 for CIE clear, 2 for CIE overcast,
+				  3 for uniform, 4 for CIE intermediate
 	A2			- zenith brightness
 	A3			- ground plane brightness
 	A4			- normalization factor based on sun direction
 	A5,A6,A7		- sun direction
 }
 
-skybright = wmean((Dz+1.01)^10, if(A1, cloudysky, sunnysky),
-		(Dz+1.01)^-10, A3);
+cosgamma = Dx*A5 + Dy*A6 + Dz*A7;
 
-wmean(a, x, b, y) = (a*x + b*y) / (a + b);
+gamma = Acos(cosgamma);		{ angle from sun to this point in sky }
 
-cloudysky = A2 * (1 + 2*Dz) / 3;
+zt = Acos(A7);			{ angle from zenith to sun }
 
+eta = Acos(Dz);			{ angle from zenith to this point in sky }
+
+wmean(a, x, b, y) : (a*x + b*y) / (a + b);
+
+skybr = wmean((Dz+1.01)^10,
+		select(A1, sunnysky, cloudysky, unifsky, intersky),
+		(Dz+1.01)^-10, A3);
+
 sunnysky = A2 * (.91 + 10*exp(-3*gamma) + .45*cosgamma*cosgamma)
 	      *  if( Dz - .01, 1.0 - exp(-.32/Dz), 1.0) / A4;
 
-gamma = Acos(cosgamma);
+cloudysky = A2 * (1 + 2*Dz)/3;
 
-cosgamma = Dx*A5 + Dy*A6 + Dz*A7;
+unifsky = A2;
+
+intersky = A2 * ( (1.35*sin(5.631-3.59*eta)+3.12)*sin(4.396-2.6*zt)
+			+ 6.37 - eta ) / 2.326 *
+		exp(gamma*-.563*((2.629-eta)*(1.562-zt)+.812)) / A4;