--- ray/src/util/disk2square.cal	2014/07/25 16:58:20	2.3
+++ ray/src/util/disk2square.cal	2024/10/30 17:25:37	2.6
@@ -1,4 +1,4 @@
-{ RCSid $Id: disk2square.cal,v 2.3 2014/07/25 16:58:20 greg Exp $ }
+{ RCSid $Id: disk2square.cal,v 2.6 2024/10/30 17:25:37 greg Exp $ }
 {
 	Convert between unit square and disk, using Shirley-Chiu mapping
 
@@ -51,19 +51,21 @@ 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_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, -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)
@@ -71,16 +73,16 @@ out_square_y = (out_square_b + 1)/2;
 	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 = -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_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;
+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),