MODULE srRotaVox;
IMPORT srBase, srMath, Math;

TYPE Voxel=srBase.Voxel;
TYPE PT=srBase.PT;
TYPE Ray=srBase.Ray;
TYPE SREAL=srBase.SREAL;

TYPE RVox*=OBJECT(Voxel);
VAR
	child: Voxel;
	phi: REAL; (* rotation angle about Z axis in radians *)
	dtick: REAL;
	center:PT;

PROCEDURE&init*(c:Voxel; dt: REAL);
BEGIN
	child:=c;
	dtick:=dt;
	register;
	center.x:=1/2; center.y:=1/2; center.z:=1/2;
END init;

PROCEDURE tick;
BEGIN
	phi:=phi+dtick;
	IF phi > 6.2832 THEN phi:=0 END
END tick;

PROCEDURE d2(a,b:PT):SREAL;
BEGIN
	 RETURN((a.x-b.x)*(a.x-b.x)+ (a.y-b.y)*(a.y-b.y) + (a.z-b.z)*(a.z-b.z));
END d2;

PROCEDURE dia(a,b:PT):SREAL;
BEGIN
	 RETURN Math.sqrt((a.x-b.x)*(a.x-b.x)+ (a.y-b.y)*(a.y-b.y) + (a.z-b.z)*(a.z-b.z));
END dia;

PROCEDURE Shade(VAR ray: Ray);
VAR
	a,b,c,dxyz:PT;
	dc,x,y,z: REAL;
	i: INTEGER;
BEGIN
	(* advance ray to its intersection with the sphere of radius 1/2 centered in the voxel*)
	a:=ray.lxyz;
	b:= srBase.Exit(ray);
	c.x := (a.x+b.x)/2; c.y := (a.y+b.y)/2; c.z := (a.z + b.z)/2;
	FOR i := 0 TO 12 DO
		dc := d2(center,c);
		IF dc < 1/4 THEN
			b:=c;
		ELSE
			a:=c;
		END;
		c.x := (a.x+b.x)/2; c.y := (a.y+b.y)/2; c.z := (a.z + b.z)/2;
	END;
	dc := d2(center,c);
	IF ABS(dc-1/4)<0.001 THEN
		a.x:=(x-1/2)*2; a.y:=(y-1/2)*2; a.z:=(z-1/2)*2;      (* we know that the vector from (1/2,1/2,1/2) to (x,y,z) has length 1/2 so this translates and normalizeds it *)
		srMath.orrot(a,Zaxis,phi);
		a.x:=a.x/2+1/2; a.y:=a.y/2+1/2; a.z:=a.z/2+1/2;	(* denormalize and translate bacck*)
		ray.lxyz:=a;
		dxyz:=ray.dxyz;
		srMath.orrot(ray.dxyz,Zaxis,-phi);
		ray.normal:=ray.dxyz;
		child.Shade(ray);
		IF ~ray.changed THEN
			ray.dxyz:=dxyz

		ELSE		(* reflection or refraction. Must rotate back to enclosing coordinate system *)
					(* 2 B implemented *)
			ray.changed := FALSE;
			ray.dxyz:=dxyz
		END
	END;
END Shade;

END RVox;

VAR
	Zaxis: PT;
BEGIN
	Zaxis.z:=1;

END srRotaVox.