--- ray/src/rt/rayinit.cal	2013/03/25 16:10:00	2.15
+++ ray/src/rt/rayinit.cal	2019/06/10 13:56:52	2.19
@@ -1,4 +1,4 @@
-{ RCSid $Id: rayinit.cal,v 2.15 2013/03/25 16:10:00 greg Exp $ }
+{ RCSid $Id: rayinit.cal,v 2.19 2019/06/10 13:56:52 greg Exp $ }
 {
 	Initialization file for Radiance.
 
@@ -37,6 +37,9 @@
 
 	select(N, a1, a2, ..)		- return aN
 
+	min(a1, a2, ..)			- return minimum argument
+	max(a1, a2, ..)			- return maximum argument
+
 	sqrt(x)				- square root function
 
 	sin(x), cos(x), tan(x),
@@ -51,8 +54,6 @@
 
 	rand(x)				- pseudo-random function (0 to 1)
 
-	hermite(p0,p1,r0,r1,t)		- 1-dimensional hermite polynomial
-
 	noise3(x,y,z), noise3x(x,y,z),
 	noise3y(x,y,z), noise3z(x,y,z)	- noise function with gradient (-1 to 1)
 
@@ -95,8 +96,6 @@ xor(a,b) : if( a, not(b), b );
 abs(x) : if( x, x, -x );
 sgn(x) : if( x, 1, if(-x, -1, 0) );
 sq(x) : x*x;
-max(a,b) : if( a-b, a, b );
-min(a,b) : if( a-b, b, a );
 inside(a,x,b) : and(x-a,b-x);
 frac(x) : x - floor(x);
 mod(n,d) : n - floor(n/d)*d;
@@ -123,6 +122,11 @@ cross(i,v1,v2) : select(i,	v1(2)*v2(3) - v1(3)*v2(2),
 fade(near_val,far_val,dist) : far_val +
 		if (16-dist, (near_val-far_val)/(1+dist*dist), 0);
 
+hermite(p0,p1,r0,r1,t) :	p0 * ((2*t-3)*t*t+1) +
+				p1 * (-2*t+3)*t*t +
+				r0 * (((t-2)*t+1)*t) +
+				r1 * ((t-1)*t*t);
+
 bezier(p1, p2, p3, p4, t) : 	p1 * (1+t*(-3+t*(3-t))) +
 				p2 * 3*t*(1+t*(-2+t)) +
 				p3 * 3*t*t*(1-t) +
@@ -147,21 +151,42 @@ turbulencez(x,y,z,s) : if( s-1.01, 0,
 
 			{ Normal distribution from uniform range (0,1) }
 
-un2`P(t) : t - (2.515517+t*(.802853+t*.010328))/
+un2`P.(t) : t - (2.515517+t*(.802853+t*.010328))/
 		(1+t*(1.432788+t*(.189269+t*.001308))) ;
-un1`P(p) : un2`P(sqrt(-2*log(p))) ;
+un1`P.(p) : un2`P.(sqrt(-2*log(p))) ;
 
-unif2norm(p) : if( .5-p, -un1`P(p), un1`P(1-p) ) ;
+unif2norm(p) : if( .5-p, -un1`P.(p), un1`P.(1-p) ) ;
 
 nrand(x) = unif2norm(rand(x));
 
 			{ Local (u,v) coordinates for planar surfaces }
-crosslen`P = Nx*Nx + Ny*Ny;
+crosslen`P. = Nx*Nx + Ny*Ny;
 			{ U is distance from projected Z-axis }
-U = if( crosslen`P - FTINY,
-		(Py*Nx - Px*Ny)/crosslen`P,
+U = if( crosslen`P. - FTINY,
+		(Py*Nx - Px*Ny)/crosslen`P.,
 		Px);
 			{ V is defined so that N = U x V }
-V = if( crosslen`P - FTINY,
-		Pz - Nz*(Px*Nx + Py*Ny)/crosslen`P,
+V = if( crosslen`P. - FTINY,
+		Pz - Nz*(Px*Nx + Py*Ny)/crosslen`P.,
 		Py);
+
+			{ Local hemisphere direction for *func & *data types }
+			{ last 3 real args = unnormalized up-vector }
+Vux`P. = arg(AC-1)*NzP - arg(AC)*NyP;
+Vuy`P. = arg(AC)*NxP - arg(AC-2)*NzP;
+Vuz`P. = arg(AC-2)*NyP - arg(AC-1)*NxP;
+vnorm`P. = 1/sqrt(Vux`P.*Vux`P. + Vuy`P.*Vuy`P. + Vuz`P.*Vuz`P.);
+Vnx`P. = Vux`P.*vnorm`P.;
+Vny`P. = Vuy`P.*vnorm`P.;
+Vnz`P. = Vuz`P.*vnorm`P.;
+Unx`P. = NyP*Vnz`P. - NzP*Vny`P.;
+Uny`P. = NzP*Vnx`P. - NxP*Vnz`P.;
+Unz`P. = NxP*Vny`P. - NyP*Vnx`P.;
+			{ Transform vectors, normalized (dx,dy,dz) away from surf }
+Ldx(dx,dy,dz) = dx*Unx`P. + dy*Uny`P. + dz*Unz`P.;
+Ldy(dx,dy,dz) = dx*Vnx`P. + dy*Vny`P. + dz*Vnz`P.;
+Ldz(dx,dy,dz) = dx*NxP + dy*NyP + dz*NzP;
+			{ Incident vector transformed to our coords }
+Idx = Ldx(-Dx,-Dy,-Dz);
+Idy = Ldy(-Dx,-Dy,-Dz);
+Idz = RdotP;