// -----------------------------------------------------------------------------
// This work is licensed under the Creative Commons Attribution-ShareAlike 3.0 Unported License.
// To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/3.0/ or send a
// letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View, California, 94041, USA.
//
// Persistence of Vision Raytracer Scene Description File
// File: Apollonian_Gasket.pov
// Description: Apollonian gasket or Apollonian net fractal
//
// 2024-05-25 - Louis bellotto - louisbel@free.fr
// -----------------------------------------------------------------------------
#version 3.8;



// -----------------------------------------------------------------------------
// --- INCLUDES ----------------------------------------------------------------
// -----------------------------------------------------------------------------
#include "colors.inc"

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



// -----------------------------------------------------------------------------
// --- GLOBAL SETTINGS ---------------------------------------------------------
// -----------------------------------------------------------------------------
global_settings {
 assumed_gamma 1.0
 }

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



// -----------------------------------------------------------------------------
// --- PARAMETERS --------------------------------------------------------------
// -----------------------------------------------------------------------------
#declare DarkMode = true;

#declare rLine = 0.02;

#declare BG_COLOR = Black;
#declare FORE_COLOR = Black;
#if(DarkMode) #declare FORE_COLOR = White; #end


#declare RecursionLevel = 10;



// -----------------------------------------------------------------------------
// --- SCENE SETTINGS ----------------------------------------------------------
// -----------------------------------------------------------------------------
camera {
 orthographic
 location <0, 0, -48>
 up y
 right x*image_width/image_height
 look_at <0, 0, 0>
 angle 36
 rotate -9*z
 }

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

#if(DarkMode)
 background { color BG_COLOR }
#else
 plane { z,0 pigment { checker White,White*0.95 } }
#end




// -----------------------------------------------------------------------------
// --- MACROS & FUNCTIONS ------------------------------------------------------
// -----------------------------------------------------------------------------

// --- FOR THE FIRST FOUR CIRCLES ----------------------------------------------
#macro DrawCircle(V,C)
 torus {
  V.z, rLine
  rotate 90*x
  pigment { color C }
  translate <V.x,V.y,-0.10>
  }
#end

// --- COMPLEX MULTIPLICATION --------------------------------------------------
// --- Use <u,v> vector of POVRay as a + ib
#macro complexMult(z1,z2)
 <z1.u*z2.u - z1.v*z2.v, z1.v*z2.u + z1.u*z2.v>
#end

// --- COMPLEX SQUARE ROOT -----------------------------------------------------
#macro complexSqrt(z1)
 #if(z1.v=0)
  #if(z1.u>0)
   #local result=<sqrt(z1.u),0>;
  #else
   #local result=<0,sqrt(-z1.u)>;
  #end
 #else
  #local sqrtr=sqrt(sqrt(z1.u*z1.u + z1.v*z1.v));
  #local theta=atan2(z1.v,z1.u); 
  #local result=<sqrtr*cos(theta/2),sqrtr*sin(theta/2)>;
 #end
 result
#end

// --- INNER SODDY CIRCLE ------------------------------------------------------
#macro SoddyCircleInternal(C1,C2,C3)
 #local a = sqrt(pow(C2.x-C3.x,2)+pow(C2.y-C3.y,2));
 #local b = sqrt(pow(C3.x-C1.x,2)+pow(C3.y-C1.y,2));
 #local c = sqrt(pow(C2.x-C1.x,2)+pow(C2.y-C1.y,2));
 #local Delta = 0.5*abs((C2.x-C1.x)*(C3.y-C1.y)-(C2.y-C1.y)*(C3.x-C1.x));

 #local tc1 = 1 + 2*Delta/(a*(b+c-a));
 #local tc2 = 1 + 2*Delta/(b*(c+a-b));
 #local tc3 = 1 + 2*Delta/(c*(a+b-c));
 #local den = a*tc1 + b*tc2 + c*tc3;

 #local k1 = a*tc1/den;
 #local k2 = b*tc2/den;
 #local k3 = c*tc3/den;

 #local r1 = C1.z; #local r2 = C2.z; #local r3=C3.z;
 < k1*C1.x + k2*C2.x + k3*C3.x, k1*C1.y + k2*C2.y + k3*C3.y, 1/(1/r1+1/r2+1/r3+2*sqrt(1/r1/r2+1/r2/r3+1/r3/r1)) >
#end

// --- OUTER SODDY CIRCLE ------------------------------------------------------
#macro SoddyCircleExternal(C1,C2,C3)
 #local r1 = C1.z;
 #local r2 = C2.z;
 #local r3=-C3.z;
 
 #local r = 1/(1/r1+1/r2+1/r3+2*sqrt(1/r1/r2+1/r2/r3+1/r3/r1));
 
 #local z1 = <C1.x,C1.y>;
 #local z2 = <C2.x,C2.y>;
 #local z3 = <C3.x,C3.y>;

 #local term1 = 1/r1*z1+1/r2*z2+1/r3*z3;
 #local term2 = 2*complexSqrt(1/r1/r2*complexMult(z1,z2)+1/r2/r3*complexMult(z2,z3)+1/r1/r3*complexMult(z1,z3));

 #local center1 = r*(term1-term2);
 #local center2 = r*(term1+term2);

 #local d1 = pow(center1.u-C3.x,2) + pow(center1.v-C3.y,2);
 #local d2 = pow(center2.u-C3.x,2) + pow(center2.v-C3.y,2);

 #if(d1>d2)
  #local out = <center1.u, center1.v, r>;
 #else
  #local out = <center2.u, center2.v, r>;
 #end

 out
#end

// --- RECURSIVE APOLLONIAN GASKET ---------------------------------------------
#macro ApollonianGasket(C1, C2, C3, n, exter)
 #if(exter)
  #local Cnew = SoddyCircleExternal(C1,C2,C3);
 #else
  #local Cnew = SoddyCircleInternal(C1,C2,C3);
 #end

 union {
  torus { Cnew.z, rLine rotate 90*x translate <Cnew.x, +Cnew.y, 0> }
  torus { Cnew.z, rLine rotate 90*x translate <Cnew.x, -Cnew.y, 0> }
  pigment { color FORE_COLOR }
  }

 #if(n>1)
  ApollonianGasket(Cnew,C2,C3,n-1,exter)
  ApollonianGasket(C1,Cnew,C3,n-1,exter)
  ApollonianGasket(C1,C2,Cnew,n-1,0)
 #end

#end



// -----------------------------------------------------------------------------
// --- STARTING CIRCLES --------------------------------------------------------
// -----------------------------------------------------------------------------
#declare C1 = < 9.750000,  0.000000,  5.250000>;
#declare C2 = <-5.250000,  0.000000,  9.750000>;
#declare C3 = < 5.825243,  8.834951,  4.417476>;
#declare C4 = < 0.000000,  0.000000, 15.000000>;

DrawCircle(C1, FORE_COLOR)
DrawCircle(C2, FORE_COLOR)
DrawCircle(C3, FORE_COLOR)
DrawCircle(C4, FORE_COLOR)




// -----------------------------------------------------------------------------
// --- RECURSIVE APOLLONIAN GASKET FRACTAL -------------------------------------
// -----------------------------------------------------------------------------
ApollonianGasket(C1, C2, C3, RecursionLevel, 0)
ApollonianGasket(C1, C2, C4, RecursionLevel, 1)
ApollonianGasket(C1, C3, C4, RecursionLevel, 1)
ApollonianGasket(C2, C3, C4, RecursionLevel, 1)




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