{ RCSid $Id: disk2square.cal,v 2.5 2019/08/13 18:29:23 greg Exp $ }
{
	Convert between unit square and disk, using Shirley-Chiu mapping

	5/20/2011	G. Ward

	Inputs:
		in_disk_x	Input unit disk Cartesian coordiantes
		in_disk_y	center at (0,0)
	or:
		in_square_x	Input unit square [0,1]^2
		in_square_y

	Corresponding outputs:
		out_square_x	Output unit square [0,1]^2
		out_square_y
	or:
		out_disk_x	Output unit disk Cartesian coordinates
		out_disk_y	center at (0,0)
}
			{ Compute disk position from square coordinates }
in_square_a = 2*in_square_x - 1;
in_square_b = 2*in_square_y - 1;

in_square_rgn = if(in_square_a + in_square_b,
			if(in_square_a - in_square_b, 1, 2),
			if(in_square_b - in_square_a, 3, 4));

out_disk_r = select(in_square_rgn, in_square_a, in_square_b,
			-in_square_a, -in_square_b);

out_disk_phi = PI/4 * select(in_square_rgn,
				in_square_b/in_square_a,
				2 - in_square_a/in_square_b,
				4 + in_square_b/in_square_a,
				if(in_square_b*in_square_b,
					6 - in_square_a/in_square_b,
					0));

out_disk_x = out_disk_r*cos(out_disk_phi);
out_disk_y = out_disk_r*sin(out_disk_phi);

			{ Compute square position from disk coordinates }
norm_radians(p) : if(-p - PI/4, p + 2*PI, p);
in_disk_r = sqrt(in_disk_x*in_disk_x + in_disk_y*in_disk_y);
in_disk_phi = norm_radians(atan2(in_disk_y, in_disk_x));

in_disk_rgn = floor((in_disk_phi + PI/4)/(PI/2)) + 1;

out_square_a = select(in_disk_rgn,
			in_disk_r,
			(PI/2 - in_disk_phi)*in_disk_r/(PI/4),
			-in_disk_r,
			(in_disk_phi - 3*PI/2)*in_disk_r/(PI/4),
			in_disk_r);

out_square_b = select(in_disk_rgn,
			in_disk_phi*in_disk_r/(PI/4),
			in_disk_r,
			(PI - in_disk_phi)*in_disk_r/(PI/4),
			-in_disk_r, -in_disk_r);

out_square_x = (out_square_a + 1)/2;
out_square_y = (out_square_b + 1)/2;
{
	The following forumulas compute Shirley-Chiu bin "scbin" based on:

	RHS		- right-handed system (-1 for left-handed coords)
	Dx,Dy,Dz	- Incident direction (normalized, towards surface front)
	rNx,rNy,rNz	- Surface normal (normalized, away from surface)
	Ux,Uy,Uz	- Up direction vector (does not need to be normalized)

	The SCdim variable assigns the square side dimension for bins, which are
	ordered with the "up" direction changing fastest.
}
RHS = 1;
				{ Compute oriented axis values }
inc_dz = -Dx*rNx-Dy*rNy-Dz*rNz;
inc_rx = -RHS*(Dx*(Uy*rNz-Uz*rNy) + Dy*(Uz*rNx-Ux*rNz) + Dz*(Ux*rNy-Uy*rNx));
inc_ry = -Dx*Ux-Dy*Uy-Dz*Uz - inc_dz*(rNx*Ux+rNy*Uy+rNz*Uz);
inc_den2 = inc_rx*inc_rx + inc_ry*inc_ry;
inc_radf = if(inc_den2-1e-7, sqrt((1 - inc_dz*inc_dz)/inc_den2), 0);
				{ Pass to formulas in first section }
in_disk_x = inc_rx*inc_radf;
in_disk_y = inc_ry*inc_radf;
				{ Compute final bin (-1 if behind surface) }
scbin = if(inc_dz,
		floor(out_square_x*SCdim)*SCdim + floor(out_square_y*SCdim),
		-1);