3D Steiner chain and Villarceau circles with POV-Ray

HierachicalSteinerChainVillarceau.md

HierachicalSteinerChainVillarceau.pov

#version 3.7;

global_settings { assumed_gamma 1 }
#default { finish{ ambient 0.1 diffuse 0.9 }}

#include "circumsphere.inc" // https://gist.github.com/stla/995699c09c49a2d6033bbd5d2c2cf628
#include "TorusThreePoints2.inc" // https://gist.github.com/stla/fffd87b9648fc6ac882ee6de75a062d3 
#include "colors.inc"
#include "textures.inc"
#include "glass.inc"
#include "metals.inc"

// Villarceau circles ----------------------------------------------------------
#macro villarceau(mu, a, c, theta, psi, epsilon)
    #local b = sqrt(a*a-c*c);
    #local bb = b*sqrt(mu*mu-c*c);
    #local bb2 = b*b*(mu*mu-c*c);
    #local denb1 = c*(a*c-mu*c+c*c-a*mu-bb);
    #local b1 = (a*mu*(c-mu)*(a+c)-bb2+c*c+bb*(c*(a-mu+c)-2*a*mu))/denb1;
    #local denb2 = c*(a*c-mu*c-c*c+a*mu+bb);
    #local b2 = (a*mu*(c+mu)*(a-c)+bb2-c*c+bb*(c*(a-mu-c)+2*a*mu))/denb2;
    #local omegaT = (b1+b2)/2;
    #local d = (a-c)*(mu-c)+bb;
    #local r = c*c*(mu-c)/((a+c)*(mu-c)+bb)/d;
    #local R = c*c*(a-c)/((a-c)*(mu+c)+bb)/d;
    #local omega2 = (a*mu + bb)/c;
    #local sign = (epsilon > 0 ? 1 : -1);
    #local f1 = -sqrt(R*R-r*r)*sin(theta);
    #local f2 = sign*(r+R*cos(theta));
    #local x1 = f1*cos(psi) + f2*sin(psi) + omegaT;
    #local y1 = f1*sin(psi) - f2*cos(psi);
    #local z1 = r*sin(theta);
    #local den = pow(x1-omega2,2)+y1*y1+z1*z1;
    < omega2 + (x1-omega2)/den, y1/den, z1/den >
#end
#declare gold = texture {
    pigment { BrightGold }
    finish {
      ambient .1
      diffuse .1
      reflection .25
      specular 1
      metallic
      //brilliance 5
    }
}
#macro vcircles(mu,a,c,n,r, shift) 
    #for(i,0,n-1)
        #local psi = 2*i*pi/n;
        #local A1 = villarceau(mu, a, c, 0, psi, 1);
        #local B1 = villarceau(mu, a, c, 2, psi, 1);
        #local C1 = villarceau(mu, a, c, 4, psi, 1);
        Torus(A1, B1, C1, r, texture{gold}, shift)
        #local A2 = villarceau(mu, a, c, 0, psi, -1);
        #local B2 = villarceau(mu, a, c, 2, psi, -1);
        #local C2 = villarceau(mu, a, c, 4, psi, -1);
        Torus(A2, B2, C2, r, texture{gold}, shift)
    #end
#end

// main macro; parameter -1 < phi < 1, phi /= 0 --------------------------------
#macro Steiner3D(n, phi, Center, Radius, Sphere, Ellipse, Cyclide, Villarceau, Depth, Sense)
    #local Sign = (Sense>0 ? 1 : -1);
    #local sine = sin(pi/n);
    #local halfSide = Radius*sine; // side length of the n-gon 
    #local Coef = 1/(1+sine); // Radius/(Radius+halfSide)
    #local invphi = 1/phi;
    #macro inversion(M,R,C)
        #local II = <invphi*R,0,0>; // inversion pole
        #local k = R*R*(invphi*invphi-1); // negated inversion constant
        #local IM = M - II - C;
        II + C - k/vdot(IM, IM)*IM
    #end
    //// ---------------- cyclide; notations Garnier & al ----------------- ////
    // center of exterior sphere 
    #local O1 = <2*invphi*Radius,0,0>;
    // interior sphere
    #if(n>2)
        #local SmallRadius = Coef*(Radius-halfSide);
        #local p1 = inversion(<SmallRadius,0,0>, Radius, 0);
        #local p2 = inversion(<0,SmallRadius,0>, Radius, 0);
        #local p3 = inversion(<-SmallRadius,0,0>, Radius, 0);
        #local p4 = inversion(<0,0,SmallRadius>, Radius, 0);
        #local cs = circumsphere(p1, p2, p3, p4);
        #local O2 = <cs.x, cs.y, cs.z>; // cs.z should be 0
        #local r = cs.t;
    #else
       #local O2 = inversion(<0,0,0>, Radius, 0);
       #local r = 0;
    #end
    // a and b are also used for the ellipse
    #local c = -(O1.x - O2.x)/2; 
    #local a = (r + Radius)/2;
    #local b = sqrt(a*a - c*c);
    #local mu = Radius - a;
    #if(Cyclide)
        polynomial {4
            xyz(2,0,0): -2*mu*mu + 2*b*b - 4*a*a,
            xyz(1,0,0): 8*a*c*mu,
            xyz(0,0,0): -4*c*c*mu*mu + mu*mu*mu*mu + b*b*b*b - 2*mu*mu*b*b,
            xyz(0,2,0): -2*mu*mu - 2*b*b,
            xyz(0,0,2): -2*mu*mu + 2*b*b,
            xyz(2,2,0): 2,
            xyz(2,0,2): 2,
            xyz(0,2,2): 2,
            xyz(4,0,0): 1,
            xyz(0,4,0): 1,
            xyz(0,0,4): 1
            texture{ Glass }
            translate Center + (O2-O1)/2 
        }
    #end
    #if(Villarceau)
        vcircles(mu, a, c, 25, 0.02, Center + (O2-O1)/2)
    #end

    // ellipse -----------------------------------------------------------------
    #if(Ellipse)
        #local nbsegs=500;
        sphere_sweep {
            linear_spline 1+nbsegs
            #for(i,0,nbsegs)
                #local ui = 2*pi*i/nbsegs;
                < a*cos(ui), b*sin(ui), 0 > 0.025
            #end
            pigment{ Red }
            translate Center + (O2-O1)/2
            no_shadow
        }
    #end

    // -------------------------------------------------------------------------
    #local shift = Sign*pi/90;
    #local ChalfSide = Coef*halfSide;
    #local CRadius = Coef*Radius;
    #local i=1;
    #while(i<=n)
        #local beta = 2*i*pi/n + frame_number*shift; // frame 1 to 180 or 180/n
        #local pti = <CRadius*cos(beta), CRadius*sin(beta), 0> + Center;
        #local p1 = inversion(<ChalfSide,0,0> + pti, Radius, Center);
        #local p2 = inversion(<0,ChalfSide,0> + pti, Radius, Center);
        #local p3 = inversion(<-ChalfSide,0,0> + pti, Radius, Center);
        #local p4 = inversion(<0,0,ChalfSide> + pti, Radius, Center);
        #local cs = circumsphere(p1, p2, p3, p4);
        #local center = <cs.x, cs.y, 0> - O1;
        #local r = cs.t;
        #if(Depth=1)
            sphere {
                center, r
                texture { Chrome_Metal }
                finish {
                    ambient .1
                 //   reflection 0.0
                    specular 1
                    roughness .1
                }
            }
        #end
        #if(Depth>1)
            Steiner3D(n, phi, center, r, 0, 0, Cyclide, 0, Depth-1, -Sense)
        #end
        #local i = i+1;
    #end

    // starting sphere ---------------------------------------------------------
    #if(Sphere)
        sphere {
            Center, Radius
            texture {
                Glass
                finish {
                    ambient 0.7
                    reflection 0.0
                    specular 1
                    roughness 0.001
                }
            }
        }
    #end
#end

// "sky" -----------------------------------------------------------------------
sky_sphere { 
    pigment { 
        gradient <0,1,0>
        color_map { 
            [0.00 rgb <0.00, 0.01, 0.05>]
            [0.35 rgb <0.00, 0.00, 0.00>]
            [0.65 rgb <0.00, 0.00, 0.00>]
            [1.00 rgb <0.00, 0.01, 0.05>]
        }
    }
}

// -----------------------------------------------------------------------------
#declare Center = <0,0,0>; // arbitrary in plane z=0
#declare Radius = 3; // arbitrary >0

// camera and light source -----------------------------------------------------
camera {
    location <2, -8, -8>
    look_at Center
    angle 40
    rotate <0,0,0>
}
light_source { <0, 0, -60> White }

// -----------------------------------------------------------------------------
#declare depth = 1;
Steiner3D(5, 0.4, Center, Radius, 0, 0, 0, 1, depth, 1)

                                                                           
/* ini file
Quality = 9
Antialias = on
Antialias_Threshold = 0.4
Antialias_Depth = 3
Antialias_Gamma = 1
Sampling_Method = 1
Jitter_Amount = 1
Bounding = on
Input_File_Name = HierarchicalSteinerChainVillarceau.pov
Height = 512
Width = 512
Initial_Frame = 1
Final_Frame = 36
Cyclic_Animation = on
*/