--- ray/src/rt/rayinit.cal	1991/08/22 08:57:19	1.1
+++ ray/src/rt/rayinit.cal	1992/10/21 11:20:26	2.5
@@ -3,16 +3,13 @@
 {
 	Initialization file for Radiance.
 
-	4/14/86
-}
-
-{
 	The following are predefined:
 
 	Dx, Dy, Dz			- ray direction
 	Nx, Ny, Nz			- surface normal
 	Px, Py, Pz			- intersection point
 	T				- distance from start
+	Ts				- single ray (shadow) distance
 	Rdot				- ray dot product
 	S				- world scale
 	Tx, Ty, Tz			- world origin
@@ -70,7 +67,7 @@ noise3d(i,x,y,z) = select(i, noise3a(x,y,z), noise3b(x
 bound(a,x,b) : if(a-x, a, if(x-b, b, x));
 Acos(x) : acos(bound(-1,x,1));
 Asin(x) : asin(bound(-1,x,1));
-Exp(x) : if(-x-60, 0, exp(x));
+Exp(x) : if(-x-100, 0, exp(x));
 Sqrt(x) : if(x, sqrt(x), 0);
 
 			{ Useful constants }
@@ -95,15 +92,15 @@ linterp(t,p0,p1) : (1-t)*p0 + t*p1;
 
 noop(v) = v;
 clip(v) = bound(0,v,1);
-noneg(v) = max(0,v);
-red(r,g,b) = r;
-green(r,g,b) = g;
-blue(r,g,b) = b;
-grey(r,g,b) = .3*r + .59*g + .11*b;
+noneg(v) = if(v,v,0);
+red(r,g,b) = if(r,r,0);
+green(r,g,b) = if(g,g,0);
+blue(r,g,b) = if(b,b,0);
+grey(r,g,b) = noneg(.263*r + .655*g + .082*b);
 clip_r(r,g,b) = bound(0,r,1);
 clip_g(r,g,b) = bound(0,g,1);
 clip_b(r,g,b) = bound(0,b,1);
-clipgrey(r,g,b) = bound(0,grey(r,g,b),1);
+clipgrey(r,g,b) = min(grey(r,g,b),1);
 
 dot(v1,v2) : v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3);
 cross(i,v1,v2) : select(i,	v1(2)*v2(3) - v1(3)*v2(2),
@@ -134,3 +131,24 @@ turbulenceb(x,y,z,s) = if( s-1.01, 0,
 turbulencec(x,y,z,s) = if( s-1.01, 0,
 			sgn(noise3(x/s,y/s,z/s))*noise3c(x/s,y/s,z/s) +
 			turbulencec(x,y,z,2*s) );
+
+			{ Normal distribution from uniform range (0,1) }
+
+un2`private(t) : t - (2.515517+t*(.802853+t*.010328))/
+		(1+t*(1.432788+t*(.189269+t*.001308))) ;
+un1`private(p) : un2`private(sqrt(-2*log(p))) ;
+
+unif2norm(p) : if( .5-p, un1`private(p), -un1`private(1-p) ) ;
+
+nrand(x) = unif2norm(rand(x));
+
+			{ Local (u,v) coordinates for planar surfaces }
+crosslen`private = Nx*Nx + Ny*Ny;
+			{ U is distance from origin in XY-plane }
+U = if( crosslen`private - FTINY,
+		(Py*Nx - Px*Ny)/crosslen`private,
+		Px);
+			{ V is defined so that N = U x V }
+V = if( crosslen`private - FTINY,
+		Pz - Nz*(Px*Nx + Py*Ny)/crosslen`private,
+		Py);