On 6/13/24 10:17, kurtz le pirate wrote:
> wip here :<http://louisbel.free.fr/scenes/scene038.shtml>

Cool!

Bill P.

Post a reply to this message

"Bald Eagle" <cre### [at] netscapenet> wrote:

> I just started watching Daniel Schiffman's video this morning - I was hoping to
> start writing some code based on it after work.

Well, I got partway there.
I'm botching _something_ up (as usual).

I need to stop blindly hacking and do some more thinking and learning.


#version 3.8;
global_settings {assumed_gamma 1.0 }

#declare S3 = sqrt(3);

#declare M = 3;
#declare C1 = < 1,  0,  1>+M*(x+y);
#declare C2 = <-1,  0,  1>+M*(x+y);
#declare C3 = < 0, S3,  1>+M*(x+y);


camera {
 location M*<1, S3/2, -2.25>
 right x*image_width/image_height
 up y
 look_at M*<1, S3/2, 0>
 //rotate y*15
}
sky_sphere {pigment {rgb 1}}
light_source {< 0, 0, -150> rgb 1}

#declare Line = 0.005;

cylinder {-x*100, x*100, Line pigment {rgb x}}
cylinder {-y*100, y*100, Line pigment {rgb y}}

#declare k4curvature1 = function (k1, k2, k3) {k1 + k2 + k3 + 2 * sqrt (k1*k2 +
k2*k3 + k3*k1)}
#declare k4curvature2 = function (k1, k2, k3) {k1 + k2 + k3 - 2 * sqrt (k1*k2 +
k2*k3 + k3*k1)}



// Make every circle a vector: <Real, Imaginary, Radius>
// Start out with <R, I, 0> to do math, and then tack on radius in last step
// That way add, subtract, and scalar multiplication are all just operations on
a vector

#macro Re(Z) Z.x #end
#macro Im(Z) Z.y #end
#macro Mult(z1, z2) <Re(z1)*Re(z2) - Im(z1)*Im(z2), Re(z1)*Im(z2) +
Im(z1)*Re(z2)> #end

#macro Sqrt (Z)
 #local m = sqrt ( Re(Z)*Re(Z) + Im(Z)*Im(Z) );
 #local Angle = atan2 (Im(Z), Re(Z)) / 2;
 #local NewZ = sqrt(m)*<cos(Angle), sin(Angle)>;
 NewZ
#end

#macro ComplexDescartes (C1, C2, C3, k4)
 // z41 = (z1*k1 + z2*k2 + z3*k3 + 2 * sqrt (z1*k1*z2*k2 + z2*k2*z3*k3 +
z3*k3*z1*k1)) / k4
 // z42 = (z1*k1 + z2*k2 + z3*k3 - 2 * sqrt (z1*k1*z2*k2 + z2*k2*z3*k3 +
z3*k3*z1*k1)) / k4
 #local k1 = 1/C1.z;
 #local k2 = 1/C2.z;
 #local k3 = 1/C3.z;

 #local z1 = <C1.x, C1.y>;
 #local z2 = <C2.x, C2.y>;
 #local z3 = <C3.x, C3.y>;

 #local zk1 = z1*k1;
 #local zk2 = z2*k2;
 #local zk3 = z3*k3;

 #local Sum = zk1 + zk2 + zk3;
 #local Square = Mult (zk1, zk2) + Mult (zk2, zk3) + Mult (zk3, zk1);
 #local Root = Sqrt (Square);

 #local z41 = (Sum + 2*Root) / k4;
 #local z42 = (Sum - 2*Root) / k4;

 #local Z41 = z41 + <0, 0, 1/k4>;
 #local Z42 = z42 + <0, 0, 1/k4>;

 #local Result = array [2] {Z41, Z42};

 Result
#end

#macro Circle (C)
 torus {abs(C.z), Line pigment {rgb 0} rotate x*90 translate <C.x, C.y, 0>}
#end








#declare k41 = k4curvature1 (1/C1.z, 1/C2.z, 1/C3.z);
#declare k42 = k4curvature2 (1/C1.z, 1/C2.z, 1/C3.z);

#declare C41 = ComplexDescartes (C1, C2, C3, k41);
#declare C42 = ComplexDescartes (C1, C2, C3, k42);

#if (0)
Circle (C1)
Circle (C2)
Circle (C3)
Circle (C41[0])
Circle (C42[1])
#end

#macro DrawCircles (Array)
 #local Circles = dimension_size (Array, 1)-1;
 #for (C, 0, Circles)
  Circle (Array [C])
 #end
#end

#declare AllCircles = array;
#declare AllCircles [0] = C1;
#declare AllCircles [1] = C2;
#declare AllCircles [2] = C3;

#declare Queue = array;
#declare Queue [0] = array [3] {C1, C2, C3};

//#declare Queue [0] = C1;
//#declare Queue [1] = C2;
//#declare Queue [2] = C3;

#declare Levels = 3;

#macro Gasket (RecursionLevels)
 #for (L, 0, RecursionLevels)
  #local Triplets = dimension_size (Queue, 1)-1;
  #for (T, 0, Triplets)
   DrawCircles (Queue [T])
  #end

  #declare NextQueue = array;
  #local QIndex = 0;
  #for (T, 0, Triplets)
   #local triplet = Queue [T];
   #local c1 = triplet [0];
   #local c2 = triplet [1];
   #local c3 = triplet [2];
   #local k41 = k4curvature1 (1/c1.z, 1/c2.z, 1/c3.z);
   #local k42 = k4curvature2 (1/c1.z, 1/c2.z, 1/c3.z);
   #declare new1 = ComplexDescartes (c1, c2, c3, k41);
   #declare new2 = ComplexDescartes (c1, c2, c3, k42);

   #local T1 = array [3] {c1, c2, new1[0]};
   #local T2 = array [3] {c2, c3, new1[1]};
   #local T3 = array [3] {c3, c1, new1[1]};
   #local T4 = array [3] {c1, c2, new2[1]};
   #local T5 = array [3] {c2, c3, new2[1]};
   #local T6 = array [3] {c3, c1, new2[1]};
   //#local T6 = array [3] {c2, c3, new1[1]};

   #declare NextQueue [QIndex] = T1;
   #local QIndex = QIndex + 1;
   #declare NextQueue [QIndex] = T2;
   #local QIndex = QIndex + 1;
   #declare NextQueue [QIndex] = T3;

   #local QIndex = QIndex + 1;
   #declare NextQueue [QIndex] = T4;
   #local QIndex = QIndex + 1;
   #declare NextQueue [QIndex] = T5;
   #local QIndex = QIndex + 1;
   #declare NextQueue [QIndex] = T6;

  #end // end for T

  #declare Queue = NextQueue;
 #end // end for L
#end

Gasket (Levels)

Post a reply to this message

Attachments:
Download 'appolonian_gasket.png' (62 KB)

Preview of image 'appolonian_gasket.png'

"Bald Eagle" <cre### [at] netscapenet> wrote:

> Well, I got partway there.
> I'm botching _something_ up (as usual).

I almost got it.
Removing one of the circle solutions cleaned up a bunch of unnecessary circles,
but the circles for the top right area are in the wrong place.
I'm wondering if there's some sign or floating point thing going on.

This is 5 levels of recursion.

Minor changes to recursion macro:
(mostly just 2 changes to an array index, restoring them to my initial values)

#macro Gasket (RecursionLevels)
 #for (L, 0, RecursionLevels)
  #local Triplets = dimension_size (Queue, 1)-1;
  #for (T, 0, Triplets)
   DrawCircles (Queue [T])
  #end

  #declare NextQueue = array;
  #local QIndex = 0;
  #for (T, 0, Triplets)
   #local triplet = Queue [T];
   #local c1 = triplet [0];
   #local c2 = triplet [1];
   #local c3 = triplet [2];
   #local k41 = k4curvature1 (1/c1.z, 1/c2.z, 1/c3.z);
   #local k42 = k4curvature2 (1/c1.z, 1/c2.z, 1/c3.z);
   #declare new1 = ComplexDescartes (c1, c2, c3, k41);
   #declare new2 = ComplexDescartes (c1, c2, c3, k42);

   #local T1 = array [3] {c1, c2, new1[0]};
   #local T2 = array [3] {c2, c3, new1[0]};//
   #local T3 = array [3] {c3, c1, new1[0]};//
   #local T4 = array [3] {c1, c2, new2[1]};
   #local T5 = array [3] {c2, c3, new2[1]};
   #local T6 = array [3] {c3, c1, new2[1]};
   //#local T6 = array [3] {c2, c3, new1[1]};

   #declare NextQueue [QIndex] = T1;
   #local QIndex = QIndex + 1;
   #declare NextQueue [QIndex] = T2;
   #local QIndex = QIndex + 1;
   #declare NextQueue [QIndex] = T3;

   #local QIndex = QIndex + 1;
   #declare NextQueue [QIndex] = T4;
   //#local QIndex = QIndex + 1;
   //#declare NextQueue [QIndex] = T5;
   #local QIndex = QIndex + 1;
   #declare NextQueue [QIndex] = T6;

  #end // end for T

  #declare Queue = NextQueue;
 #end // end for L
#end

Post a reply to this message

Attachments:
Download 'appolonian_gasket.png' (58 KB)

Preview of image 'appolonian_gasket.png'

"Bald Eagle" <cre### [at] netscapenet> wrote:

> I recall in one of the sites that I stumbled across, there was a method for
> calculating how many circles there were for any given level of recursion.
> That might be a nice addition to the debug stream, as well as warning/sanity
> check before a render is attempted to be started with too high a level of
> recursion for the system it's being run on.

https://en.wikibooks.org/wiki/Fractals/Apollonian_fractals

Post a reply to this message

A few quick explanations for those interested.
(I'd probably do a page on my site... but I'm so far behind ;) )


I haven't used Bald Eagle's inversion method, which looks very
interesting and promising.


The reference being the WIKIPEDIA page
<https://en.wikipedia.org/wiki/Pappus_chain>



-- The two basic radii :
#declare BigRadius = SIZE; // AB/2
#declare SmallRadius = SIZE*3/4; // AC/2


-- For the first circle (CB):
#declare cx = SmallRadius*2 + BigRadius - SmallRadius;
#declare r = BigRadius - SmallRadius;
DrawCircle(cx, 0, r)


-- Next, a simple loop :
#declare r = SmallRadius/BigRadius;
#declare NthRadius = r;

#declare n = 1;
#while ( NthRadius > MIN_RADIUS )
  #declare divisor = n*n*(1-r)*(1-r)+r;
  #declare xx = 2*BigRadius*r*(1+r)/(2*divisor);
  #declare yy = 2*BigRadius*n*r*(1-r)/divisor;
  #declare NthRadius = 2*BigRadius*r*(1-r)/(2*divisor);

  DrawCircle(xx, +yy, NthRadius)
  DrawCircle(xx, -yy, NthRadius) // symmetrically

  #declare n = n + 1;
#end


The first tests were in 2D. So I made a macro that drew circles :
#macro DrawCircle(X,Y,R,C)
 torus { R, rLine pigment { color C } rotate 90*x translate <X,Y,0> }
#end


Transformed for 3d :
#macro DrawCircle(X,Y,R,C)
 sphere { <X,R,Y>, R pigment { color C } }
#end



A yes, I forgot! The end of the loop is controlled by the constancy
MIN_RADIUS. If the radius is "too" small, we stop.

#declare MIN_RADIUS = 0.50; for the stone spheres





-- 
Kurtz le pirate
Compagnie de la Banquise

Post a reply to this message

On 14/06/2024 16:30, Bald Eagle wrote:
> "Bald Eagle" <cre### [at] netscapenet> wrote:
> https://en.wikibooks.org/wiki/Fractals/Apollonian_fractals

This guy is working on similar things and making interactive toys:
https://observablehq.com/@esperanc/3d-apollonian-sphere-packings
--
YB

Post a reply to this message

On 14/06/2024 19:30, yesbird wrote:
> This guy is working on similar things and making interactive toys:
> https://observablehq.com/@esperanc/3d-apollonian-sphere-packings

Also this paper may be of interest:
https://arxiv.org/pdf/1609.03811
--
YB

Post a reply to this message

On 12/06/2024 17:52, Bald Eagle wrote:
> I'll bet that would amazing done in glass with caustics.
> Or as stone spheres in a field of grass.
> 

After grass and stones, glass.




-- 
Kurtz le pirate
Compagnie de la Banquise

Post a reply to this message

Attachments:
Download 'papus4_colors_1.png' (349 KB)

Preview of image 'papus4_colors_1.png'

Ugh.

I tried rewriting from scratch (porting the Javascript code from
https://en.wikibooks.org/wiki/Fractals/Apollonian_fractals to SDL)
and I'm still stuck with that upper left hand portion not being right.

I thought maybe some of the quadratic solution checking in that could would fix
things, because saving that code as an SVG file and opening it in a browser
works beautifully (and it's FAST!).

So, still sitting here with my dunce cap on until I figure out how to wrap my
head around what's going wrong with the center coordinates and the recursion
queue.

- BW

Post a reply to this message

Am 16.06.2024 um 13:44 schrieb Bald Eagle:
> Ugh.
> 
> I tried rewriting from scratch (porting the Javascript code from
> https://en.wikibooks.org/wiki/Fractals/Apollonian_fractals to SDL)
> and I'm still stuck with that upper left hand portion not being right.
> 
> I thought maybe some of the quadratic solution checking in that could would fix
> things, because saving that code as an SVG file and opening it in a browser
> works beautifully (and it's FAST!).
> 
> So, still sitting here with my dunce cap on until I figure out how to wrap my
> head around what's going wrong with the center coordinates and the recursion
> queue.
> 
> - BW
> 
It took me a while to understand this weird javascript recursion and I 
mixed up some - with + or vice versa:

global_settings {assumed_gamma 1.0 }

#include "colors.inc"

#declare S3 = sqrt(3);

#declare M = 3;
#declare C1 = < 1,  0,  1>+M*(x+y);
#declare C2 = <-1,  0,  1>+M*(x+y);
#declare C3 = < 0, S3,  1>+M*(x+y);


camera {
  location M*<0,0,-2.25>//M*<1, S3/2, -2.25>
  right x*image_width/image_height
  up y
  look_at <0,0,0>//M*<1, S3/2, 0>
  //rotate y*15
}
sky_sphere {pigment {rgb 1}}
light_source {< 0, 0, -150> rgb 1}

#declare Line = 0.005;

cylinder {-x*100, x*100, Line pigment {rgb x}}
cylinder {-y*100, y*100, Line pigment {rgb y}}

#declare k4curvature1 = function (k1, k2, k3) {k1 + k2 + k3 + 2 * sqrt 
(k1*k2 +
k2*k3 + k3*k1)}
#declare k4curvature2 = function (k1, k2, k3) {k1 + k2 + k3 - 2 * sqrt 
(k1*k2 +
k2*k3 + k3*k1)}

#macro Re(Z) Z.x #end
#macro Im(Z) Z.y #end
#macro Cmul(z1, z2)
    <Re(z1)*Re(z2) - Im(z1)*Im(z2), Re(z1)*Im(z2) + Im(z1)*Re(z2)>
#end
#macro Cadd(a,b)
    <Re(a)+Re(b),Im(a)+Im(b)>
#end
#macro Csub(a,b)
    <Re(a)-Re(b),Im(a)-Im(b)>
#end

#macro Csqrt (Z)
  #local m = sqrt ( Re(Z)*Re(Z) + Im(Z)*Im(Z) );
  #local Angle = atan2 (Im(Z), Re(Z)) / 2;
  #local NewZ = sqrt(m)*<cos(Angle), sin(Angle)>;
  NewZ
#end
/*#macro Csqrt(Z)
    #local m = sqrt ( Re(Z)*Re(Z) + Im(Z)*Im(Z) );
    <sqrt((m+Re(Z))/2),abs(Im(Z))/Im(Z)*sqrt((m-Re(Z))/2)>
#end*/


#declare Colours=array[7]{Violet,Blue,Cyan,Green,Yellow,Orange,Red}

#macro DrawCircle (C,PigmNr)
    torus {abs(C.z), Line pigment {colour Colours[PigmNr]} rotate x*90 
translate <C.x, C.y, 0>}
#end


/***********************************************
Edit the seed for other circles or run an animation
#declare Zufall=seed(36730+frame_number);
************************************************/

#declare Zufall=seed(36730+17);

// unit circle
#declare b=<0,0,-1>;

// insert two raqndomly positioned touching circles
#declare tr = 1-rand(Zufall)/2;
#declare pa = pi/2 - asin(rand(Zufall)*(1-tr)/tr);
#declare px = tr * cos(pa);
#declare py = tr * sin(pa);
#declare pr = 1 - tr;
#declare qr = (1 - pr) * (1 - cos(pa + pi/2))  / (1 + pr - (1 - pr) 
*cos(pa + pi/2));
#declare qx = 0;
#declare qy = qr - 1;

#declare p=<px,py,pr>;
#declare q=<qx,qy,qr>;

#declare MinRadius=0.1;

#declare MaxDepth=7;

#macro Kiss(a,b,c,initial,Depth)
    #if (Depth <= MaxDepth)
       #local k1 = 1/a.z; // real
       #local z1 = <a.x,a.y>; // complex
       #local kz1 = Cmul(<k1,0>,z1); // complex

       #local k2 = 1/b.z; // real
       #local z2 = <b.x,b.y>; // complex
       #local kz2 = Cmul(<k2,0>,z2); // complex

       #local k3 = 1/c.z; // real
       #local z3 = <c.x,c.y>; // complex
       #local kz3 = Cmul(<k3,0>,z3); // complex

       #local k4p = k1 + k2 + k3 + 2*sqrt(k1*k2 + k2*k3 + k3*k1); // real
       #local k4m = k1 + k2 + k3 - 2*sqrt(k1*k2 + k2*k3 + k3*k1); // real

       #local kz4p = 
Cadd(Cadd(Cadd(kz1,kz2),kz3),Cmul(<2,0>,Csqrt(Cadd(Cadd(Cmul(kz1,kz2),Cmul(kz2,kz3)),Cmul(kz3,kz1)))));

// complex
       #local kz4m = 
Csub(Cadd(Cadd(kz1,kz2),kz3),Cmul(<2,0>,Csqrt(Cadd(Cadd(Cmul(kz1,kz2),Cmul(kz2,kz3)),Cmul(kz3,kz1)))));

// complex

       #if (k4p > k4m)
          #local k4 = k4p; // real
          #local kz4 = kz4p; // complex
          #local k4b = k4m; // real
          #local kz4b = kz4m; // complex
       #else
          #local k4 = k4m; // real
          #local kz4 = kz4m; // complex
          #local k4b = k4p; // real
          #local kz4b = kz4p; // complex
       #end

       #local cc = <kz4.x/k4 , kz4.y/k4 , abs(1/k4)>; // Circle

       #local dx = a.x - cc.x;
       #local dy = a.y - cc.y;
       #local dr = a.z + cc.z;

       #local dtest=abs(dx*dx + dy*dy - dr*dr);

       #if ( abs(dx*dx + dy*dy - dr*dr) > 0.0001 )
          #local cc = <kz4b.x/k4 , kz4b.y/k4 , abs(1/k4)>;
       #end

       #if (initial)
          #local cc=<kz4b.x/k4b,kz4b.y/k4b,abs(1/k4b)>;
       #end

       DrawCircle(cc,mod(Depth,7))

       Kiss(a,b,cc,false,Depth+1)
       Kiss(a,cc,c,false,Depth+1)
       Kiss(cc,b,c,false,Depth+1)
    #end
#end


DrawCircle(b,0)
DrawCircle(p,0)
DrawCircle(q,0)


Kiss(b,p,q,true,1)
Kiss(b,p,q,false,1)

Post a reply to this message

Attachments:
Download 'apollonian_fractals_03.png' (68 KB)

Preview of image 'apollonian_fractals_03.png'