POV-Ray : Newsgroups : povray.advanced-users : Pappus Chain : Re: Pappus Chain Server Time18 Jul 2024 08:26:21 EDT (-0400)
 Re: Pappus Chain
 From: MichaelJF Date: 17 Jun 2024 15:53:19 Message: <6670942f@news.povray.org>
```
{
"@context": "https://schema.org",
"@type": "DiscussionForumPosting",
"@id": "#6670942f%40news.povray.org",
"dateCreated": "2024-06-17T19:53:19+00:00",
"datePublished": "2024-06-17T19:53:19+00:00",
"author": {
"@type": "Person",
"name": "MichaelJF"
}
}
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
<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 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 =

// complex
#local kz4m =

// 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)
```

Attachments: