We assume that the radius of the outer stationary circle is **1**. If we set the radius **x** of the inner stationary circle, all the other circles are uniquely determined by solving the system **Sys**. Should be **x<=1/3** . If **x=1/3** then all the inner circles have a radius **1/3** . The following picture explains the meaning of symbols in the procedure **Circles**:

**Circles:=proc(x)**

**local OO, O1, O2, O3, O4, O2x, O2y, O3x, O3y, OT, T1, T2, T3, s, t, dist, Sys, Sol, sol, y, u, v, z, C0, R0, P;**

**uses plottools, plots;**

**OO:=[0,0]: O1:=[x+y,0]: O2:=[O2x,O2y]: O3:=[O3x,O3y]: O4:=[-x-z,0]: OT:=[x+2*y-1,0]:**

**T1:=(O2*~y+O1*~u)/~(y+u): T2:=(O3*~u+O2*~v)/~(u+v): T3:=(O4*~v+O3*~z)/~(v+z):**

**solve({(T2-T1)[1]*(s-((T1+T2)/2)[1])+(T2-T1)[2]*(t-((T1+T2)/2)[2])=0, (T3-T2)[1]*(s-((T2+T3)/2)[1])+(T3-T2)[2]*(t-((T3+T2)/2)[2])=0}, {s,t}):**

**assign(%);**

**dist:=(A,B)->sqrt((B[1]-A[1])^2+(B[2]-A[2])^2):**

**Sys:={dist(O1,O2)^2=(y+u)^2, dist(OO,O2)^2=(x+u)^2, dist(O2,O3)^2=(u+v)^2, dist(OO,O3)^2=(x+v)^2, dist(O3,O4)^2=(z+v)^2, x+y+z=1, dist(O2,OT)^2=(1-u)^2, dist(O3,OT)^2=(1-v)^2};**

**Sol:=op~([allvalues([solve(Sys)])]);**

**sol:=select(i->is(eval(convert([y>0,u>0,v>0,z>0,O2y>0,x<=y,u<=y,v<=u,z<=v],`and`),i)), Sol)[];**

**assign(sol);**

**O1:=[x+y,0]: O2:=[O2x,O2y]: O3:=[O3x,O3y]: O4:=[-x-z,0]: OT:=[x+2*y-1,0]:**

**C0:=eval([s,t],sol);**

**R0:=eval(dist(T1,C0),sol):**

**P:=proc(phi)**

**local eq, r1, r, R, Ot, El, i, S, s, t, P1, P2;**

**uses plots,plottools;**

**eq:=1-dist([r*cos(s),r*sin(s)],OT)=r-x;**

**r1:=solve(eq,r);**

**r:=eval(r1,s=phi);**

**R[1]:=evalf(r-x);**

**Ot[1]:=evalf([r*cos(phi),r*sin(phi)]);**

**El:=plot([r1*cos(s),r1*sin(s),s=0..2*Pi],color="Green",thickness=3);**

**for i from 2 to 6 do**

**S:=[solve({1-dist(OT,[s,t])=dist(Ot[i-1],[s,t])-R[i-1], 1-dist(OT,[s,t])=dist(OO,[s,t])-x})];**

**P1:=eval([s,t],S[1]); P2:=eval([s,t],S[2]);**

**Ot[i]:=`if`(evalf(Ot[i-1][1]*P1[2]-Ot[i-1][2]*P1[1])>0,P1,P2);**

**R[i]:=dist(Ot[i],OO)-x;**

**od;**

**display(El,seq(disk(Ot[k],0.012),k=1..6),circle(C0,R0,color=gold,thickness=3),circle([x+2*y-1,0],1, color=blue,thickness=4), circle(OO,x, color=red,thickness=4), seq(circle(Ot[k],R[k], thickness=3),k=1..6), scaling=constrained, axes=none);**

**end proc:**

**animate(P,[phi], phi=0..Pi, frames=120);**

**end proc: **

Example of use (I got **x=0.22** just by measuring the ruler displayed original animation):

**Circles(0.22);**

The curve on the following animation is **an astroid **(a special case of **hypocycloid**). See wiki for details. **Hypocycloid** procedure creates animation for any hypocycloid. Parameters of the procedure: **R** is the radius of the outer circle, **r** is the radius of the inner circle.

**Hypocycloid:=proc(R,r)**

**local A, B, f, g, F;**

**uses plots,plottools;**

**A:=circle(R,color=green,thickness=4):**

**B:=display(circle([R-r,0],r,color=red,thickness=4),line([R-r,0],[R,0],color=red,thickness=4)):**

**f:=t->plot([(R-r)*cos(s)+r*cos((R-r)/r*s),(R-r)*sin(s)-r*sin((R-r)/r*s),s=0..t],color=blue,thickness=4):**

**g:=t->rotate(rotate(B,-R/r*t,[R-r,0]),t):**

**F:=t->display(A,f(t),g(t),scaling=constrained):**

**animate(F,[t], t=0..2*Pi*denom(R/r), frames=90);**

**end proc:**

Examples of use:

**Hypocycloid(4,1); **

** **

**Hypocycloid(5,3);**

** **

Круги.mw