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'