#version 3.7;
#include "math.inc"
#include "colors.inc"

global_settings
{
	assumed_gamma	1.0
	adc_bailout		0.005
	max_trace_level	50
	ambient_light	1/3
/*
	radiosity
	{
		pretrace_start	0.08
		pretrace_end	0.01
		count		50
		error_bound	0.1
		recursion_limit	1
		normal		on
		brightness	0.8
		always_sample	yes
		gray_threshold	0.8
		media		on
	}
*/
}

camera
{
	#local cam_distance	= 80;	// render at 8192x8192
	#local cam_planesize	= 8;	// render at 8192x8192
//	#local cam_distance	= 20;	// render at 4096x4096
//	#local cam_planesize	= 2;	// render at 4096x4096
	#local cam_aspectratio	= image_width/image_height;
//	orthographic
	location	-z * cam_distance
	direction	+z * cam_distance
	right		+x * cam_planesize*cam_aspectratio
	up		+y * cam_planesize
	rotate		+x * asind(tand(30))
	rotate		-y * 045
	rotate		+y * 180
}

light_source
{
	<-32,+32,-32,>
	color rgb	1
	rotate		y * 330
	parallel
	point_at	0
	shadowless
}

light_source
{
	<-32,+32,-32,>
	color rgb	1
	rotate		y * 090
	parallel
	point_at	0
//	shadowless
}

background {srgb 1/2}


#declare XYZWhiteReference1 =  95.047;
#declare XYZWhiteReference2 = 100.000;
#declare XYZWhiteReference3 = 108.883;
#declare XYZEpsilon = 0.008856;
#declare XYZKappa = 903.3;


// input L = between 0 and 100
// input C = between 0 and 100
// input H = between 0 and 360
// output L = between 0 and 100
// output A = between -128 and +128
// output B = between -128 and +128
#declare funcLCH2LABa1 = function(L,C,H) {L}
#declare funcLCH2LABa2 = function(L,C,H) {cos(radians(H)) * C}
#declare funcLCH2LABa3 = function(L,C,H) {sin(radians(H)) * C}


// input L = between 0 and 100
// input A = between -128 and +128
// input B = between -128 and +128
// output X: between 0 and 100
// output Y: between 0 and 100
// output Z: between 0 and 100
// Note that some of these functions are out of order.
#declare funcLAB2XYZa2 = function(L,A,B) {(L + 16) / 116}
#declare funcLAB2XYZa1 = function(L,A,B) {funcLAB2XYZa2(L,A,B) + A / 500}
#declare funcLAB2XYZa3 = function(L,A,B) {funcLAB2XYZa2(L,A,B) - B / 200}
#declare funcLAB2XYZb1 = function(L,A,B) {XYZWhiteReference1 * select(pow(funcLAB2XYZa1(L,A,B), 3) - XYZEpsilon, (funcLAB2XYZa1(L,A,B) - 16.0 / 116.0) / 7.787, (funcLAB2XYZa1(L,A,B) - 16.0 / 116.0) / 7.787, pow(funcLAB2XYZa1(L,A,B), 3))}
#declare funcLAB2XYZb2 = function(L,A,B) {XYZWhiteReference2 * select(L - XYZKappa * XYZEpsilon, L / XYZKappa, L / XYZKappa, pow(((L + 16) / 116), 3))}
#declare funcLAB2XYZb3 = function(L,A,B) {XYZWhiteReference3 * select(pow(funcLAB2XYZa3(L,A,B), 3) - XYZEpsilon, (funcLAB2XYZa3(L,A,B) - 16.0 / 116.0) / 7.787, (funcLAB2XYZa3(L,A,B) - 16.0 / 116.0) / 7.787, pow(funcLAB2XYZa3(L,A,B), 3))}


// input X: between 0 and 100
// input Y: between 0 and 100
// input Z: between 0 and 100
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
// Note that out-of-range colors are *NOT* corrected. You will have to do this yourself.
#declare funcXYZ2RGBa1 = function(X,Y,Z) {X/100 *  3.2406 + Y/100 * -1.5372 + Z/100 * -0.4986}
#declare funcXYZ2RGBa2 = function(X,Y,Z) {X/100 * -0.9689 + Y/100 *  1.8758 + Z/100 *  0.0415}
#declare funcXYZ2RGBa3 = function(X,Y,Z) {X/100 *  0.0557 + Y/100 * -0.2040 + Z/100 *  1.0570}
#declare funcXYZ2RGBb1 = function(X,Y,Z) {select(funcXYZ2RGBa1(X,Y,Z) - 0.0031308, 12.92 * funcXYZ2RGBa1(X,Y,Z), 12.92 * funcXYZ2RGBa1(X,Y,Z), 1.055 * pow(funcXYZ2RGBa1(X,Y,Z), 1 / 2.4) - 0.055)}
#declare funcXYZ2RGBb2 = function(X,Y,Z) {select(funcXYZ2RGBa2(X,Y,Z) - 0.0031308, 12.92 * funcXYZ2RGBa2(X,Y,Z), 12.92 * funcXYZ2RGBa2(X,Y,Z), 1.055 * pow(funcXYZ2RGBa2(X,Y,Z), 1 / 2.4) - 0.055)}
#declare funcXYZ2RGBb3 = function(X,Y,Z) {select(funcXYZ2RGBa3(X,Y,Z) - 0.0031308, 12.92 * funcXYZ2RGBa3(X,Y,Z), 12.92 * funcXYZ2RGBa3(X,Y,Z), 1.055 * pow(funcXYZ2RGBa3(X,Y,Z), 1 / 2.4) - 0.055)}


// input R: between 0 and 1
// input G: between 0 and 1
// input B: between 0 and 1
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
#declare correctRGB1 = function(R,G,B) {min(max(R,0),1)}
#declare correctRGB2 = function(R,G,B) {min(max(G,0),1)}
#declare correctRGB3 = function(R,G,B) {min(max(B,0),1)}


#local VLCH = <50, 50, 180>;
#local VLAB = <funcLCH2LABa1(VLCH.x,VLCH.y,VLCH.z), funcLCH2LABa2(VLCH.x,VLCH.y,VLCH.z), funcLCH2LABa3(VLCH.x,VLCH.y,VLCH.z)>;
#local VXYZ = <funcLAB2XYZb1(VLAB.x,VLAB.y,VLAB.z), funcLAB2XYZb2(VLAB.x,VLAB.y,VLAB.z), funcLAB2XYZb3(VLAB.x,VLAB.y,VLAB.z)>;
#local VRGB = <funcXYZ2RGBb1(VXYZ.x,VXYZ.y,VXYZ.z), funcXYZ2RGBb2(VXYZ.x,VXYZ.y,VXYZ.z), funcXYZ2RGBb3(VXYZ.x,VXYZ.y,VXYZ.z)>;
#local VRGB = <correctRGB1(VRGB.x,VRGB.y,VRGB.z), correctRGB2(VRGB.x,VRGB.y,VRGB.z), correctRGB3(VRGB.x,VRGB.y,VRGB.z)>;
#debug concat("\nVLCH = ", vstr(3, VLCH, ", ", 0, 5), "\n")
#debug concat("\nVLAB = ", vstr(3, VLAB, ", ", 0, 5), "\n")
#debug concat("\nVXYZ = ", vstr(3, VXYZ, ", ", 0, 5), "\n")
#debug concat("\nVRGB = ", vstr(3, VRGB, ", ", 0, 5), "\n")

#local VLCH1 = 50;
#local VLCH2 = 50;
#local VLCH3 = 180;
#local VLAB1 = funcLCH2LABa1(VLCH1,VLCH2,VLCH3);
#local VLAB2 = funcLCH2LABa2(VLCH1,VLCH2,VLCH3);
#local VLAB3 = funcLCH2LABa3(VLCH1,VLCH2,VLCH3);
#local VXYZ1 = funcLAB2XYZb1(VLAB1,VLAB2,VLAB3);
#local VXYZ2 = funcLAB2XYZb2(VLAB1,VLAB2,VLAB3);
#local VXYZ3 = funcLAB2XYZb3(VLAB1,VLAB2,VLAB3);
#local VRGB1 = funcXYZ2RGBb1(VXYZ1,VXYZ2,VXYZ3);
#local VRGB2 = funcXYZ2RGBb2(VXYZ1,VXYZ2,VXYZ3);
#local VRGB3 = funcXYZ2RGBb3(VXYZ1,VXYZ2,VXYZ3);
#local VRGB1 = correctRGB1(VRGB1,VRGB2,VRGB3);
#local VRGB2 = correctRGB2(VRGB1,VRGB2,VRGB3);
#local VRGB3 = correctRGB3(VRGB1,VRGB2,VRGB3);

#declare convertLCH2RGB_R = function(L,C,H) {funcXYZ2RGBb1(funcLAB2XYZb1(funcLCH2LABa1(L,C,H),funcLCH2LABa2(L,C,H),funcLCH2LABa3(L,C,H)),funcLAB2XYZb2(funcLCH2LABa1(L,C,H),funcLCH2LABa2(L,C,H),funcLCH2LABa3(L,C,H)),funcLAB2XYZb3(funcLCH2LABa1(L,C,H),funcLCH2LABa2(L,C,H),funcLCH2LABa3(L,C,H)))}
#declare convertLCH2RGB_G = function(L,C,H) {funcXYZ2RGBb2(funcLAB2XYZb1(funcLCH2LABa1(L,C,H),funcLCH2LABa2(L,C,H),funcLCH2LABa3(L,C,H)),funcLAB2XYZb2(funcLCH2LABa1(L,C,H),funcLCH2LABa2(L,C,H),funcLCH2LABa3(L,C,H)),funcLAB2XYZb3(funcLCH2LABa1(L,C,H),funcLCH2LABa2(L,C,H),funcLCH2LABa3(L,C,H)))}
#declare convertLCH2RGB_B = function(L,C,H) {funcXYZ2RGBb3(funcLAB2XYZb1(funcLCH2LABa1(L,C,H),funcLCH2LABa2(L,C,H),funcLCH2LABa3(L,C,H)),funcLAB2XYZb2(funcLCH2LABa1(L,C,H),funcLCH2LABa2(L,C,H),funcLCH2LABa3(L,C,H)),funcLAB2XYZb3(funcLCH2LABa1(L,C,H),funcLCH2LABa2(L,C,H),funcLCH2LABa3(L,C,H)))}


isosurface
{
	function {convertLCH2RGB_R(y*200-100,sqrt(x*x+z*z)*100,atan2(x,z)*pi/180)}
	contained_by
	{
		box {-2,+2}
	}
	threshold 1
	accuracy	0.01
	max_gradient	5
//	[evaluate P0, P1, P2]
//	open
//	[max_trace INTEGER] | [all_intersections]
	pigment {rgb 1}
}