// Title: LCH Color Solid Cylinder Isosurface
// Author: Michael Horvath, Christoph Lipka
// Website: http://isometricland.net
// Created: 2016-11-20
// Updated: 2016-11-28
// This file is licensed under the terms of the CC-GNU LGPL.
// http://www.gnu.org/licenses/lgpl-2.1.html
// Illuminant =  D65
// Observer   =  2° (1931)
// +kfi0 +kff15 +kc

#version 3.7;
#include "math.inc"
#include "colors.inc"
#declare ScaleAmount = 200;


global_settings
{
	assumed_gamma	1.0
	adc_bailout	0.005
	max_trace_level	50
	ambient_light	1/3
/*
	ambient_light	0
	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	= 40;	// render at 8192x8192
//	#local cam_planesize	= 4;	// render at 8192x8192
	#local cam_distance	= 21;
	#local cam_planesize	= 2.1;
//	#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		+x * 90
	rotate		-y * 045
//	rotate		+y * 180
	scale		ScaleAmount
}

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

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

//background {srgb 1}


#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 128
// 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) {cosd(H) * C}
#declare funcLCH2LABa3 = function(L,C,H) {sind(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)}


// utility
#declare selectRGBa = function(C,D) {select(C,-1,D)}
#declare selectRGBb = function(C,D) {select(C-1,D,D,-1)}
#declare linearRGBa = function(C) {C/12.92}
#declare linearRGBb = function(C) {pow((C + 0.055)/(1 + 0.055), 2.4)}
#declare linearRGBc = function(C) {select(C-0.04045,linearRGBa(C),linearRGBa(C),linearRGBb(C))}


// 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
// only uses these when converting individual values, not inside isosurfaces or function patterns
#declare correctRGBa1 = function(R,G,B) {selectRGBa(R,selectRGBa(G,selectRGBa(B,R)))}
#declare correctRGBa2 = function(R,G,B) {selectRGBa(R,selectRGBa(G,selectRGBa(B,G)))}
#declare correctRGBa3 = function(R,G,B) {selectRGBa(R,selectRGBa(G,selectRGBa(B,B)))}
#declare correctRGBb1 = function(R,G,B) {selectRGBb(R,selectRGBb(G,selectRGBb(B,correctRGBa1(R,G,B))))}
#declare correctRGBb2 = function(R,G,B) {selectRGBb(R,selectRGBb(G,selectRGBb(B,correctRGBa2(R,G,B))))}
#declare correctRGBb3 = function(R,G,B) {selectRGBb(R,selectRGBb(G,selectRGBb(B,correctRGBa3(R,G,B))))}


// input L = between 0 and 100
// input C = between 0 and 128
// input H = between 0 and 360
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
#declare convertLCH2RGBa1 = 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 convertLCH2RGBa2 = 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 convertLCH2RGBa3 = 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)))}
#declare convertLCH2RGBb1 = function(L,C,H) {correctRGBb1(convertLCH2RGBa1(L,C,H),convertLCH2RGBa2(L,C,H),convertLCH2RGBa3(L,C,H))}
#declare convertLCH2RGBb2 = function(L,C,H) {correctRGBb2(convertLCH2RGBa1(L,C,H),convertLCH2RGBa2(L,C,H),convertLCH2RGBa3(L,C,H))}
#declare convertLCH2RGBb3 = function(L,C,H) {correctRGBb3(convertLCH2RGBa1(L,C,H),convertLCH2RGBa2(L,C,H),convertLCH2RGBa3(L,C,H))}


// testing
#local VLCH = <25, 50, 270>;
#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 VCOR = < correctRGBb1(VRGB.x,VRGB.y,VRGB.z),  correctRGBb2(VRGB.x,VRGB.y,VRGB.z),  correctRGBb3(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")
#debug concat("\nVCOR = ", vstr(3, VCOR, ", ", 0, 5), "\n")

//#local VCON = <convertLCH2RGBa1(50, 50, 90),convertLCH2RGBa2(50, 50, 90),convertLCH2RGBa3(50, 50, 90)>;
//#local VCON = <convertLCH2RGBb1(50, 50, 180),convertLCH2RGBb2(50, 50, 180),convertLCH2RGBb3(50, 50, 180)>;
//#debug concat("\nVCON = ", vstr(3, VCON, ", ", 0, 5), "\n")


#declare fD = function(C) {abs(C-0.5)-0.5}
#declare fDist = function(Dr,Dg,Db) {max(Dr,Dg,Db)}
#declare fR = function(L,C,H) {convertLCH2RGBa1(L,C,H)}
#declare fG = function(L,C,H) {convertLCH2RGBa2(L,C,H)}
#declare fB = function(L,C,H) {convertLCH2RGBa3(L,C,H)}
//#declare fL = function(x,y,z) {acosd(y/sqrt(x*x+y*y+z*z))/180*100}
#declare fL = function(x,y,z) {atan2d(sqrt(x*x+z*z),y)/180*100}
#declare fC = function(x,y,z) {sqrt(x*x+y*y+z*z)*128}
#declare fH = function(x,y,z) {atan2d(x,z)}

#declare fFinal = function(L,C,H)
{
	fDist
	(
		fD(fR(L,C,H)),
		fD(fG(L,C,H)),
		fD(fB(L,C,H))
	)
}

#declare pigmentR = pigment
{
	function {linearRGBc(fR(fL(x,y,z),fC(x,y,z),fH(x,y,z)))}
	color_map
	{
		[0 color rgb <0,0,0>]
		[1 color rgb <3,0,0>]
	}
}
#declare pigmentG = pigment
{
	function {linearRGBc(fG(fL(x,y,z),fC(x,y,z),fH(x,y,z)))}
	color_map
	{
		[0 color rgb <0,0,0>]
		[1 color rgb <0,3,0>]
	}
}
#declare pigmentB = pigment
{
	function {linearRGBc(fB(fL(x,y,z),fC(x,y,z),fH(x,y,z)))}
	color_map
	{
		[0 color rgb <0,0,0>]
		[1 color rgb <0,0,3>]
	}
}
#declare pigmentRGB = pigment
{
	average
	pigment_map
	{
		[1 pigmentR]
		[1 pigmentG]
		[1 pigmentB]
	}
}
#declare pigmentCyl = pigment
{
	cylindrical
	color_map
	{
		[ 0/20 color rgbt 10/10]
		[ 1/20 color rgbt 10/10]
		[ 1/20 color rgbt  9/10]
		[20/20 color rgbt  9/10]
	}
	frequency 5
	scale 1.0001
}
#declare pigmentGradY = pigment
{
	gradient y
	pigment_map
	{
		[ 0/20 pigmentCyl]
		[ 1/20 pigmentCyl]
		[ 1/20 color rgbt  9/10]
		[20/20 color rgbt  9/10]
	}
	frequency 5
	scale 1.0001
}


#declare isoShape = isosurface
{
	function {fFinal(fL(x,y,z),fC(x,y,z),fH(x,y,z))}
	threshold	0
	accuracy	0.001
	contained_by
	{
		box {<-1,0,-1>,<+1,+1,+1>}
	}
	max_gradient	50	// was 20000
}

#declare sphShape = sphere
{
	0, 1
//	hollow
}

object
{
	isoShape
	pigment {pigmentRGB}
	finish {emission 0.3 diffuse 1}
	scale ScaleAmount
	rotate +y * 360 * clock
}

object
{
	sphShape
//	pigment {pigmentGradY}
	pigment {color rgbt  7/8}
	scale ScaleAmount
}