// -----------------------------------------------------------------------------
// This work is licensed under the Creative Commons Attribution 3.0 Unported License.
// To view a copy of this license, visit http://creativecommons.org/licenses/by/3.0/
// or send a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
// California, 94041, USA.
//
// Persistence Of Vision raytracer sample file.
// Demo scene : Solving Apollonius's problem
//
// 2024-06-18 - louis bellotto - louisbel@free.fr
// -----------------------------------------------------------------------------
#version 3.8;




// -----------------------------------------------------------------------------
// --- INCLUDES ----------------------------------------------------------------
// -----------------------------------------------------------------------------
#include "colors"
#include "math"




// -----------------------------------------------------------------------------
// --- VARIABLES ---------------------------------------------------------------
// -----------------------------------------------------------------------------
#if (image_width != image_height)
 #error "*** Set image ratio to 1 ! ***"
#end




// -----------------------------------------------------------------------------
#declare displayAxis = false;
#declare rLine = 0.025*2;




// -----------------------------------------------------------------------------
// --- SETTINGS ----------------------------------------------------------------
// -----------------------------------------------------------------------------
global_settings {
 assumed_gamma 1.0
 max_trace_level 5
 }

#default {
 finish {
  ambient 0.00
  emission 0.40
  diffuse 0.60
  }
 } 




// -----------------------------------------------------------------------------
// --- SCENE -------------------------------------------------------------------
// -----------------------------------------------------------------------------
#declare CameraDistance = 60;
#declare CameraAngle = 36;
#declare SIZE = int(CameraDistance*tand(CameraAngle*0.50));

camera {
 orthographic
 location <0, 0, -CameraDistance>
 up y
 right x*image_width/image_height
 look_at <0, 0, 0>
 angle CameraAngle
 }

light_source {
 <0, 0, -500>
 color White
 }

background {
 color White*0
 }

#if (displayAxis)
 #declare ArrowLen = 0.50;
 #declare FontName="times-it";
 #declare FontScale = 1.0;
 #declare Col = White;
 union {
  // x-axis
  cylinder { <-SIZE,0,0>, <+SIZE,0,0>, 0.025 }
  cylinder { <SIZE,0,0>, <SIZE-ArrowLen,+ArrowLen*tand(20),0> 0.050 }
  cylinder { <SIZE,0,0>, <SIZE-ArrowLen,-ArrowLen*tand(20),0> 0.050 }
  #declare xStr = text { ttf FontName "x" 0.25 <0, 0, 0> scale FontScale }
  object { xStr translate <SIZE-1.0, 0.30, -0.5> }
  // y-axis
  cylinder { <0, -SIZE, 0>, <0, +SIZE, 0>, 0.025 }
  cylinder { <0, +SIZE, 0>, <-ArrowLen*tand(20), SIZE-ArrowLen, 0>, 0.050 }
  cylinder { <0, +SIZE, 0>, <+ArrowLen*tand(20), SIZE-ArrowLen, 0>, 0.050 }
  #declare yStr = text { ttf FontName "y" 0.25 <0, 0, 0> scale FontScale }
  object { yStr translate <0.30, SIZE-0.7, -0.5> }
 
  pigment { color Col }
  }
#end




// -----------------------------------------------------------------------------
// --- FUNCTIONS ---------------------------------------------------------------
// -----------------------------------------------------------------------------
#declare fnline = function (a,b,x) { a*x+b }




// -----------------------------------------------------------------------------
// --- MACROS ------------------------------------------------------------------
// -----------------------------------------------------------------------------
#macro DrawCircle (Circle, Color, k, DrawCenter)
 torus {
  Circle.z, rLine*k
  rotate 90*x
  pigment { color Color }
  translate <Circle.x, Circle.y,0>
  }
 #if(DrawCenter)
  union {
   cylinder { <0,0,0>, <0,0,-0.10>, rLine*5 }
   cylinder { <-0.5, 0, 0>,<+0.5, 0, 0>, rLine }
   cylinder { <0, -0.5, 0>,<0, +0.5, 0>, rLine }
   translate <Circle.x, Circle.y,0>
   pigment { color Color }
   } 
 #end
#end




// -----------------------------------------------------------------------------
#macro NameCircle(Circle, Name, Color, Offset, FontName, FontScale)
 #local Str = text { ttf FontName Name 0.10 <0, 0, 0> scale FontScale };
 #local StrMax = max_extent(Str);
 #local Str = object { Str translate -StrMax*0.50 }
 object { Str pigment { color Color } translate <Circle.x, Circle.y,0>+Offset }
#end




// -----------------------------------------------------------------------------
#macro DrawLine(p1, p2, Color)
 cylinder { p1 ,p2, rLine pigment { color Color } }
#end




// -----------------------------------------------------------------------------
#macro LineParams(p1,p2)
 #local a = (p1.y - p2.y)/(p1.x - p2.x);
 #local b = p1.y - a*p1.x;
 (a,b)
#end




// -----------------------------------------------------------------------------
// --- OBJECTS -----------------------------------------------------------------
// -----------------------------------------------------------------------------
// circle object : <xpos, ypos, radius>

#declare Color = White*0.10;

#declare c1 = <10,8,8>;
DrawCircle(c1, Color, 1, true)
NameCircle(c1, "c1", Color, <1, 1, 0>, "times-it", 1.0)

#declare c2 = <-8,4,6>;
DrawCircle(c2, Color, 1, true)
NameCircle(c2, "c2", Color, <-1, 1, 0>, "times-it", 1.0)

#declare c3 = <6,-8,4>;
DrawCircle(c3, Color, 1, true)
NameCircle(c3, "c3", Color, <0, -1, 0>, "times-it", 1.0)




// -----------------------------------------------------------------------------
/*
#declare Color1 = White*0.75;
#declare m = 18;

// --- l1 : line c1/c2
#declare (l1a,l1b) = LineParams(c1,c2);
#declare l1x1 = m;
#declare l1y1 = fnline(l1a, l1b, l1x1);
#declare l1x2 = -m;
#declare l1y2 = fnline(l1a, l1b, l1x2);
DrawLine(<l1x1, l1y1, 0>, <l1x2, l1y2, 0>, Color1) 

// --- l2 : line c1/c3
#declare (l2a,l2b) = LineParams(c1,c3);
#declare l2x1 = m-5;
#declare l2y1 = fnline(l2a, l2b, l2x1);
#declare l2x2 = m-14;
#declare l2y2 = fnline(l2a, l2b, l2x2);
DrawLine(<l2x1, l2y1, 0>, <l2x2, l2y2, 0>, Color1) 

// --- l3 : line c2/c3
#declare (l3a,l3b) = LineParams(c2,c3);
#declare l3x1 = m-4;
#declare l3y1 = fnline(l3a, l3b, l3x1);
#declare l3x2 = -m;
#declare l3y2 = fnline(l3a, l3b, l3x2);
DrawLine(<l3x1, l3y1, 0>, <l3x2, l3y2, 0>, Color1) 
*/




// -----------------------------------------------------------------------------
// https://en.wikipedia.org/wiki/Problem_of_Apollonius#Algebraic_solutions
// -----------------------------------------------------------------------------
#macro ApolloniusSolver(c1, s1, c2, s2, c3, s3)

 #local p1 = (2*c2.y-2*c1.y)/(2*c2.x-2*c1.x);
 #local p2 = (c1.x*c1.x-c2.x*c2.x+c1.y*c1.y-c2.y*c2.y-c1.z*c1.z+c2.z*c2.z)/(2*c2.x-2*c1.x);
 #local p3 = (2*s2*c2.z-2*s1*c1.z)/(2*c2.x-2*c1.x);

 #local q1 = (2*c3.y-2*c2.y)/(2*c3.x-2*c2.x)-p1;
 #local q2 = (c2.x*c2.x-c3.x*c3.x+c2.y*c2.y-c3.y*c3.y-c2.z*c2.z+c3.z*c3.z)/(2*c3.x-2*c2.x)-p2;
 #local q3 = (2*s3*c3.z-2*s2*c2.z)/(2*c3.x-2*c2.x)-p3;
 
 #local p = -q2/q1; #local q = q3/q1;
 #local m = -p1*p-p2; #local n = p3-p1*q;

 // quadratic sol
 #local a = n*n+q*q-1;
 #local b = 2*m*n-2*n*c1.x+2*p*q-2*q*c1.y+2*s1*c1.z;
 #local c = c1.x*c1.x+m*m-2*m*c1.x+p*p+c1.y*c1.y-2*p*c1.y-c1.z*c1.z;
 #local d = b*b-4*a*c;
 
 #local rad = (-b-sqrt(d))/(2*a);
 <m+n*rad, p+q*rad, rad> 
#end




// -----------------------------------------------------------------------------
// https://www.rapidtables.com/convert/color/hsv-to-rgb.html
// -----------------------------------------------------------------------------
#macro hsv2rgb(H, S, V)
 #local H = (H<0?0:H);
 #local H = (H>=360?359:H);

 #local C = V * S;
 #local X = C*(1-abs(mod((H/60),2)-1));
 #local m = V - C;
  
 #if ((H >= 0)&(H < 60))
  #local R = C;
  #local G = X;
  #local B = 0;
 #elseif((H >= 60)&(H < 120))
  #local R = X;
  #local G = C;
  #local B = 0;
 #elseif((H >= 120)&(H < 180))
  #local R = 0;
  #local G = C;
  #local B = X;
 #elseif((H >= 180)&(H < 240))
  #local R = 0;
  #local G = X;
  #local B = C;
 #elseif((H >= 240)&(H < 300))
  #local R = X;
  #local G = 0;
  #local B = C;
 #elseif((H >= 300)&(H < 360))
  #local R = C;
  #local G = 0;
  #local B = X;
 #else
  #error "hsv2rgb convertion error"
 #end
 <(R+m),(G+m),(B+m)>
#end




// -----------------------------------------------------------------------------
#declare N = 8;
#declare Signs = array [N][3] {
 // internal : +1, external : -1
 {-1, -1, -1},
 {-1, -1, +1},
 {-1, +1, -1},
 {-1, +1, +1},
 {+1, -1, -1},
 {+1, -1, +1},
 {+1, +1, -1},
 {+1, +1, +1}
 }


// -----------------------------------------------------------------------------
#declare ColorStep = 360/(N+1);
#debug "          x        y         r\n"
#declare i = 0;
#while ( i < N)
 #declare c = ApolloniusSolver(c1, Signs[i][0], c2,  Signs[i][1], c3,  Signs[i][2]);
 #debug concat(str(i+1,2,0)," : ",str(c.x,8,4),", ",str(c.y,8,4),", ",str(c.z,8,4),"\n")
 #declare Color = hsv2rgb((i+1)*ColorStep, 1.0, 1.0);
 DrawCircle(c, Color, 0.50, true)
 NameCircle(c, concat("s",str(i+1,0,0)), Color, <0, -0.85, -0.5>, "times-it", 1.0)
 #declare i = i+1;
#end




// --- eof ---------------------------------------------------------------------