Hi Ed,
Thx for answer.
....
>>
>>"If the center of the second circle is outside the first, then the -
>>sign corresponds to externally tangent circles and the + sign to
>>internally tangent circles. " from :
>>http://mathworld.wolfram.com/TangentCircles.html
>>
>>Do I have misunderstood it or is it a bug ?
>
> It appears to be an error on the wolfram page. .....
Thx. I have send them this information.
>>eq13: (x1-x3)^2 + (y1 -y3)^2 - (rad1+rad3)^2; eq23: (x2-x3)^2 + (y2
>>-y3)^2 - (rad2+rad3)^2;
>>
>>eq23: expand(eq23);
>>eq13: expand(eq13);
>
> You may be doing it for some other reason, but Maxima doesn't require
> these two expand() statements to solve these equations.
OK
Below is the last version of Maxima code.
( I have also Lisp version )
It works, but center c3p3 ( step 3, peripheral gap nr 3) is not good.
Adam
/* ====================== code ========================
computes and draws to png file apollonian gasket
initial configuration : (k0,k1,k1,k1)
- three equal mutually tangent circles (curvature =k0)
- one outer circles ( curvature = k1)
with help of Ed Beroset
*/
/* ---------------------------- definitions
-----------------------------------------------*/
/* gives a list which is used by plot2d for drawing circle */
plot_circle(c,r):=[parametric, realpart(c)+r*cos(t), imagpart(c)+r*sin
(t)];
/* circle Descates theorem ( for curvatures ). It can also be used for
complex Descates theorem ( for centers) */
solve_n_eq(k_1,k_2,k_3):=float(ratsimp((k_1+k_2+k_3 - 2*sqrt(k_1*k_2 +
k_2*k_3 + k_1*k_3))));/* outer circle */
solve_p_eq(k_1,k_2,k_3):=float(ratsimp((k_1+k_2+k_3 + 2*sqrt(k_1*k_2 +
k_2*k_3 + k_1*k_3))));/* inner circle */
/* -----------------------------*/
rad1:100;
k1:1/rad1;
x1a:-rad1;
y1a:rad1;
c1a:x1a+y1a*%i;
x1b:-x1a;
y1b:y1a;
c1b:x1b+y1b*%i;
/* find c1c ( third circle mutually tangent to c1a and c1b ) */
eq13: (x1a-x1c)^2 + (y1a -y1c)^2 - (rad1+rad1)^2;
eq23: (x1b-x1c)^2 + (y1b -y1c)^2 - (rad1+rad1)^2;
b:solve([eq23,eq13],[x1c,y1c]);
c1c:float(rhs(b[1][1]) +rhs(b[1][2])*%i);
/* outer circle */
k0:solve_n_eq(k1,k1,k1);
rad0:abs(1/k0); /* negative curvature */
c0:solve_n_eq(k1*c1a,k1*c1b,k1*c1c)/k0;
/* inner soddy circle */
/* gap in the middle = between 3 equal circles */
k3a:solve_p_eq(k1,k1,k1);
rad3a:abs(1/k3a);
c3a:solve_p_eq(k1*c1a,k1*c1b,k1*c1c)/k3a;
/* 3 peripheral gaps between 2 inner circles and outer circle */
k3p:solve_p_eq(k1,k1,k0);
rad3p:abs(1/k3p);
c3p1:solve_p_eq(k1*c1a,k1*c1b,k0*c0)/k3p;
c3p2:solve_p_eq(k1*c1b,k1*c1c,k0*c0)/k3p;
c3p3:solve_p_eq(k1*c1a,k1*c1c,k0*c0)/k3p;
/* --------------------------- drawing code -----------------------*/
plot_list:[];
/* ------------------------------- step 0 = outer circle --------- */
plot_list:cons(plot_circle(c0,rad0),plot_list);
/* ----------------- step 1 = three equal mutually tangent circles */
plot_list:cons(plot_circle(c1a,rad1),plot_list);
plot_list:cons(plot_circle(c1b,rad1),plot_list);
plot_list:cons(plot_circle(c1c,rad1),plot_list);
/* --------------------------- - step 2 -------------------------*/
plot_list:cons(plot_circle(c3a,rad3a),plot_list); /* gap in the middle */
/* 3 peripheral gaps between 2 inner circles and outer circle */
plot_list:cons(plot_circle(c3p1,rad3p),plot_list);
plot_list:cons(plot_circle(c3p2,rad3p),plot_list);
plot_list:cons(plot_circle(c3p3,rad3p),plot_list);
set_plot_option([plot_format, gnuplot]);
set_plot_option([nticks, 80]);
set_plot_option([t,-%pi,%pi]);
set_plot_option([gnuplot_preamble, "unset key"])$ /* */
plot2d (plot_list)$
/*========================== end ==================================*/