Apollonius圆和反复计算得到的图案

标签：mathematica   维基   

这种问题首选的地址是维基百科： http://en.wikipedia.org/wiki/Circles_of_Apollonius

从作图的角度则要看笛卡尔定理（也是维基），尤其是其复数版本，同时给出半径和圆心位置的计算公式：http://en.wikipedia.org/wiki/Descartes%27_theorem

有了计算公式，就可以从任意一组（三个）Apollonius圆反复计算它们的公切圆，以及新的三圆组的公切圆，得到意想不到的图案：

1.首先，从Steiner chain (去维基！http://en.wikipedia.org/wiki/Steiner_chain）中的三个圆出发时得到的图案：

2. 从某个一个大圆和它直径上两边的两个半径1/2的小圆出发得到的结果：

圆里面数字对应直径或半径的比例；所以，可以放缩其数字（鉴于都是六的倍数）；

用Mathematica完成这些图像制作非常容易：

d[A_, B_] := Sqrt[(A - B).(A - B)];
tangentCircle[{Circle[{x1_, y1_}, r1_], Circle[{x2_, y2_}, r2_], 
    Circle[{x3_, y3_}, r3_]}] := 
  Module[{sols, x, y, r}, 
    sols = NSolve[{(x - x1)^2 + (y - y1)^2 == (r - 
           r1)^2, (x - x2)^2 + (y - y2)^2 == (r + 
           r2)^2, (x - x3)^2 + (y - y3)^2 == (r + r3)^2}, {x, y, r}];
    Circle[{x, y}, r] /. sols] /; 
   r1 == r2 + r3 && r1 > d[{x1, y1}, {x2, y2}];

tangentCircle[{Circle[{x1_, y1_}, r1_], Circle[{x2_, y2_}, r2_], 
    Circle[{x3_, y3_}, r3_]}] := 
  Module[{a2, b2, c2, d2, a3, b3, c3, d3, x, y, r, sign}, 
   If[r1 > d[{x1, y1}, {x2, y2}], sign = -1, sign = 1];
   a2 = 2 (x1 - x2); b2 = 2 (y1 - y2); c2 = 2 (sign r1 - r2);
   d2 = (x1^2 + y1^2 - r1^2) - (x2^2 + y2^2 - r2^2);
   a3 = 2 (x1 - x3); b3 = 2 (y1 - y3); c3 = 2 (sign r1 - r3);
   d3 = (x1^2 + y1^2 - r1^2) - (x3^2 + y3^2 - r3^2);
   x = (b3 d2 - b2 d3 - b3 c2 r + b2 c3 r)/(a2 b3 - b2 a3) // N;
   y = (-a3 d2 + a2 d3 + a3 c2 r - a2 c3 r)/(a2 b3 - a3 b2) // N;
   r = Min[
     Abs[r /. Solve[(x - x1)^2 + (y - y1)^2 == (r + sign r1)^2, r]]];
   Circle[{x, y}, r]];
tangentCircles[{Circle[{x1_, y1_}, r1_], Circle[{x2_, y2_}, r2_], 
    Circle[{x3_, y3_}, r3_]}] := 
  Module[{a2, b2, c2, d2, a3, b3, c3, d3, x, y, r, sign}, 
   If[r1 > d[{x1, y1}, {x2, y2}], sign = -1, sign = 1];
   a2 = 2 (x1 - x2); b2 = 2 (y1 - y2); c2 = 2 (sign r1 - r2);
   d2 = (x1^2 + y1^2 - r1^2) - (x2^2 + y2^2 - r2^2);
   a3 = 2 (x1 - x3); b3 = 2 (y1 - y3); c3 = 2 (sign r1 - r3);
   d3 = (x1^2 + y1^2 - r1^2) - (x3^2 + y3^2 - r3^2);
   x = (b3 d2 - b2 d3 - b3 c2 r + b2 c3 r)/(a2 b3 - b2 a3) // N;
   y = (-a3 d2 + a2 d3 + a3 c2 r - a2 c3 r)/(a2 b3 - a3 b2) // N;
   Circle[{x, y}, Abs[r]] /. 
    Solve[(x - x1)^2 + (y - y1)^2 == (r + sign r1)^2, r]];
triplets[{a_, b_, c_}, d_] := {{a, b, d}, {a, c, d}, {b, c, d}};
triplets[{a_, b_, c_Circle}, l_List] := triplets[{a, b, c}, #] & /@ l;
apollonianStep[{a_Circle, b_, c_}] := 
  triplets[{a, b, c}, tangentCircle[{a, b, c}]];
apollonianStep[l_List] := apollonianStep /@ l;

(*{circ1,circ2,circ3,circ4}={Circle[{0,0},1/6],Circle[{1/6-1/11,0},1/11],Circle[{-8/105,-2/35},1/14],Circle[{-3/50,2/25},1/15]};
init={{circ1,circ2,circ3},{circ1,circ2,circ4},{circ1,circ3,circ4},{circ2,circ3,circ4}};*)

{circ1, circ2, circ3} = {Circle[{0, 0}, 1/6], 
   Circle[{-1/12, 0}, 1/12], Circle[{1/12, 0}, 1/12]};
init = {{circ1, circ2, circ3}};

circs = TimeConstrained[
   Union[Flatten[
     init //. {a_, b_, Circle[p_, r_]} :> 
       apollonianStep[{a, b, Circle[p, r]}] /; r > .01]], 10];
moreCircs = 
  Union[Flatten[
    init //. {a_, b_, Circle[p_, r_]} :> 
      apollonianStep[{a, b, Circle[p, r]}] /; r > .0005]];

number[Circle[c_, r_]] := 
  Text[Style[Round[1/r/6], FontSize -> 550 r], c];
Graphics[{Blue, circs, moreCircs, 
  Map[number, DeleteCases[circs, Circle[{0, 0}, _]]]}, 
 ImageSize -> Full]

代码绝大部分来自于：

Mathematica.stackexchange.com: http://mathematica.stackexchange.com/questions/11016/is-there-a-way-to-solve-the-apollonius-circle-problem-in-mathematica

Apollonius圆和反复计算得到的图案

标签：mathematica   维基   

原文地址：http://blog.csdn.net/stereohomology/article/details/45742879