> with(plots):

> epi:=plot([4*cos(t)-cos(4*t),4*sin(t)-sin(4*t),t=0..2*Pi],color=black):
delto:=plot([2*cos(t)+cos(2*t),2*sin(t)-sin(2*t),t=0..2*Pi],color=red):

> x:=t->2*cos(t)+cos(2*t):
y:=t->2*sin(t)-sin(2*t):

xp:=D(x):yp:=D(y):
d:=12:
L:=NULL:
for n to 20 do
alpha:=n*Pi/10:
M:=NULL:
for i to 40 do
t:=i*Pi/20;
# { Les coefficients a1 et b1 dirigent la normale au rayon réfléchi. }
a1:=sin(t+alpha);b1:=cos(t+alpha);
a2:=1*b1;b2:=1*(-a1);
c1:=x(t)*a1+y(t)*b1;
c2:=(xp(t)-y(t))*(-b2)+(x(t)+yp(t))*a2;
# { Partie universelle pour les enveloppes de droites. }
delta:=a1*b2-a2*b1;
xa:=(c1*b2-c2*b1)/delta;ya:=(a1*c2-a2*c1)/delta;
M:=M,plot([[xa,ya],[x(t),y(t)],[x(t)+d*cos(alpha),y(t)+d*sin(alpha)]],color=blue);
od:
L:=L,display([M,plot([[-cos(2*alpha),sin(2*alpha)]],style=point,color=black,symbol=circle)])
od:

> display([plot([cos(u),sin(u),u=0..2*Pi],color=green),delto,epi,display([L],insequence=true)],axes=none,scaling=constrained,view=[-6..6,-6..6]);

> display([delto,L[3]],axes=none,scaling=constrained,view=[-6..6,-6..6]);