{ SCCSid "$SunId$ LBL" }

{
	Initialization file for Radiance.

	The following are predefined:

	Dx, Dy, Dz			- ray direction
	Nx, Ny, Nz			- surface normal
	Px, Py, Pz			- intersection point
	T				- distance from start
	Rdot				- ray dot product
	S				- world scale
	Tx, Ty, Tz			- world origin
	Ix, Iy, Iz			- world i unit vector
	Jx, Jy, Jz			- world j unit vector
	Kx, Ky, Kz			- world k unit vector
	arg(n)				- real arguments, arg(0) is count

	For brdf functions, the following are also available:

 	NxP, NyP, NzP			- perturbed surface normal
 	RdotP				- perturbed ray dot product
 	CrP, CgP, CbP			- perturbed material color

	Library functions:

	if(a, b, c)			- if a positive, return b, else c

	select(N, a1, a2, ..)		- return aN

	sqrt(x)				- square root function

	sin(x), cos(x), tan(x),
	asin(x), acos(x),
	atan(x), atan2(y,x)		- standard trig functions

	floor(x), ceil(x)		- g.l.b. & l.u.b.

	exp(x), log(x), log10(x)	- exponent and log functions

	erf(z), erfc(z)			- error functions

	rand(x)				- pseudo-random function (0 to 1)

	hermite(p0,p1,r0,r1,t)		- 1-dimensional hermite polynomial

	noise3(x,y,z), noise3a(x,y,z),
	noise3b(x,y,z), noise3c(x,y,z)	- noise function with gradient (-1 to 1)

	fnoise3(x,y,z)			- fractal noise function (-1 to 1)
}

			{ Backward compatibility }
AC = arg(0);
A1 = arg(1); A2 = arg(2); A3 = arg(3); A4 = arg(4); A5 = arg(5);
A6 = arg(6); A7 = arg(7); A8 = arg(8); A9 = arg(9); A10 = arg(10);

			{ Forward compatibility (?) }
D(i) = select(i, Dx, Dy, Dz);
N(i) = select(i, Nx, Ny, Nz);
P(i) = select(i, Px, Py, Pz);
noise3d(i,x,y,z) = select(i, noise3a(x,y,z), noise3b(x,y,z), noise3c(x,y,z));

			{ More robust versions of library functions }
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));
Sqrt(x) : if(x, sqrt(x), 0);

			{ Useful constants }
PI : 3.14159265358979323846;
DEGREE : PI/180;
FTINY : 1e-7;

			{ Useful functions }
and(a,b) : if( a, b, a );
or(a,b) : if( a, a, b );
not(a) : if( a, -1, 1 );
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;
tri(n,d) : abs( d - mod(n-d,2*d) );
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;
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);

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),
				v1(3)*v2(1) - v1(1)*v2(3),
				v1(1)*v2(2) - v1(2)*v2(1));

fade(near_val,far_val,dist) = far_val +
		if (16-dist, (near_val-far_val)/(1+dist*dist), 0);

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) +
				p4 * t*t*t ;

bspline(pp, p0, p1, pn, t) =	pp * (1/6+t*(-.5+t*(.5-1/6*t))) +
				p0 * (2/3+t*t*(-1+.5*t)) +
				p1 * (1/6+t*(.5+t*(.5-.5*t))) +
				pn * (1/6*t*t*t) ;

turbulence(x,y,z,s) = if( s-1.01, 0, abs(noise3(x/s,y/s,z/s)*s) +
						turbulence(x,y,z,2*s) );
turbulencea(x,y,z,s) = if( s-1.01, 0,
			sgn(noise3(x/s,y/s,z/s))*noise3a(x/s,y/s,z/s) +
			turbulencea(x,y,z,2*s) );
turbulenceb(x,y,z,s) = if( s-1.01, 0,
			sgn(noise3(x/s,y/s,z/s))*noise3b(x/s,y/s,z/s) +
			turbulenceb(x,y,z,2*s) );
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(log(1/p/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);