> restart:

> with(plots):

1. PAVAGE DE PENROSE

> R:=[seq(exp(I*Pi*j/5),j=0..9)];
phi:=(1+sqrt(5))/2;

R := [1, exp(1/5*I*Pi), exp(2/5*I*Pi), exp(3/5*I*Pi), exp(4/5*I*Pi), -1, exp(-4/5*I*Pi), exp(-3/5*I*Pi), exp(-2/5*I*Pi), exp(-1/5*I*Pi)]

phi := 1/2+1/2*5^(1/2)

On repère les tuiles de Penrose par des triplets [a,b,c], avec :
        (i) a valant 1 lorsque la tuile est un cerf-volant, et 2 lorsqu'il s'agit d'une flèche,
        (ii) b est le complexe correspondant au sommet de la tuile dont l'angle est le plus obtus,
        (iii) c est le complexe correspondant au sommet opposé à b.
La procédure transforme effectue la transformation de Penrose :

> transforme:=proc(t)
if t[1]=1

then return([1,(phi*t[2]+t[3])/(1+phi),t[3]+(t[2]-t[3])*R[10]],[1,(phi*t[2]+t[3])/(1+phi),t[3]+(t[2]-t[3])*R[2]],[2,t[3]+R[10]/(1+phi)*(t[2]-t[3]),t[3]]);

else return([1,t[2],t[3]],[2,t[3]+R[2]*(t[2]-t[3]),t[3]+R[2]*phi*(t[2]-t[3])]);

end if;

end proc;

transforme := proc (t) if t[1] = 1 then return [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[10]], [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[2]], [2, t[3]+R[10]*(t[2]-t[3])/(1+phi), t[3]] else...transforme := proc (t) if t[1] = 1 then return [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[10]], [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[2]], [2, t[3]+R[10]*(t[2]-t[3])/(1+phi), t[3]] else...transforme := proc (t) if t[1] = 1 then return [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[10]], [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[2]], [2, t[3]+R[10]*(t[2]-t[3])/(1+phi), t[3]] else...transforme := proc (t) if t[1] = 1 then return [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[10]], [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[2]], [2, t[3]+R[10]*(t[2]-t[3])/(1+phi), t[3]] else...transforme := proc (t) if t[1] = 1 then return [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[10]], [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[2]], [2, t[3]+R[10]*(t[2]-t[3])/(1+phi), t[3]] else...transforme := proc (t) if t[1] = 1 then return [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[10]], [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[2]], [2, t[3]+R[10]*(t[2]-t[3])/(1+phi), t[3]] else...transforme := proc (t) if t[1] = 1 then return [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[10]], [1, (phi*t[2]+t[3])/(1+phi), t[3]+(t[2]-t[3])*R[2]], [2, t[3]+R[10]*(t[2]-t[3])/(1+phi), t[3]] else...

La procédure penrose effectue n transformations de Penrose (la configuration initiale est constituée de 5 cerfs-volants placés circulairement autour d'un point) :

> penrose:=proc(n)
local L,M,i,t;

L:=[seq([1,R[2*i-1],0],i=1..5)];

for i from 1 to n do

        M:=L;L:=[];

        while M<>[] do

                L:=[op(L),transforme(M[1])];M:=[op(2..nops(M),M)];

        end do;

end do;

return L;

end proc;

penrose := proc (n) local L, M, i, t; L := [seq([1, R[2*i-1], 0], i = 1 .. 5)]; for i to n do M := L; L := []; while M <> [] do L := [op(L), transforme(M[1])]; M := [op(2 .. nops(M), M)] end do end do...penrose := proc (n) local L, M, i, t; L := [seq([1, R[2*i-1], 0], i = 1 .. 5)]; for i to n do M := L; L := []; while M <> [] do L := [op(L), transforme(M[1])]; M := [op(2 .. nops(M), M)] end do end do...penrose := proc (n) local L, M, i, t; L := [seq([1, R[2*i-1], 0], i = 1 .. 5)]; for i to n do M := L; L := []; while M <> [] do L := [op(L), transforme(M[1])]; M := [op(2 .. nops(M), M)] end do end do...penrose := proc (n) local L, M, i, t; L := [seq([1, R[2*i-1], 0], i = 1 .. 5)]; for i to n do M := L; L := []; while M <> [] do L := [op(L), transforme(M[1])]; M := [op(2 .. nops(M), M)] end do end do...penrose := proc (n) local L, M, i, t; L := [seq([1, R[2*i-1], 0], i = 1 .. 5)]; for i to n do M := L; L := []; while M <> [] do L := [op(L), transforme(M[1])]; M := [op(2 .. nops(M), M)] end do end do...penrose := proc (n) local L, M, i, t; L := [seq([1, R[2*i-1], 0], i = 1 .. 5)]; for i to n do M := L; L := []; while M <> [] do L := [op(L), transforme(M[1])]; M := [op(2 .. nops(M), M)] end do end do...

Il ne reste plus qu'à tracer le résultat :

> dessin:=proc(L)
local D,i;

D:=[];

for i from 1 to nops(L) do

        if L[i][1]=1

        then D:=[op(D),POLYGONS([[evalf(Re(L[i][3])),evalf(Im(L[i][3]))],[evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))),evalf(Im(L[i][3]+R[10]*(L[i][2]-L[i][3])))],[evalf(Re(L[i][2])),evalf(Im(L[i][2]))],[evalf(Re(L[i][3]+R[2]*(L[i][2]-L[i][3]))),evalf(Im(L[i][3]+R[2]*(L[i][2]-L[i][3])))]])];

        else D:=[op(D),POLYGONS([[evalf(Re(L[i][3])),evalf(Im(L[i][3]))],[evalf(Re(L[i][3]+R[10]*phi*(L[i][2]-L[i][3]))),evalf(Im(L[i][3]+R[10]*phi*(L[i][2]-L[i][3])))],[evalf(Re(L[i][2])),evalf(Im(L[i][2]))],[evalf(Re(L[i][3]+R[2]*phi*(L[i][2]-L[i][3]))),evalf(Im(L[i][3]+R[2]*phi*(L[i][2]-L[i][3])))]])];

        end if;

end do;

PLOT(op(D),AXESSTYLE(NONE),SCALING(CONSTRAINED));

end proc;

dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...dessin := proc (L) local D, i; D := []; for i to nops(L) do if L[i][1] = 1 then D := [op(D), POLYGONS([[evalf(Re(L[i][3])), evalf(Im(L[i][3]))], [evalf(Re(L[i][3]+R[10]*(L[i][2]-L[i][3]))), evalf(Im(L...

> dessin(penrose(0));

[Plot]

> dessin(penrose(1));

[Plot]

> dessin(penrose(2));

[Plot]

> dessin(penrose(3));

[Plot]

> dessin(penrose(4));

[Plot]

> dessin(penrose(6));

[Plot]

2. ALGORITHME DE TRACÉ DE PAVAGES DU PLAN

On stocke les informations contenues dans l'arborescence fondamentale dans un petit tableau. Chaque classe d'équivalence de groupe cristallographique est représentée par une ligne contenant :
        (i) son nom,
        (ii) une information sur les possibilités du réseau des translations : si la deuxième entrée est 0 (resp. 1, resp. 2), il s'agit de n'importe quel réseau (resp. du réseau hexagonal, resp. du réseau carré),
        (iii) la transformation permettant de remonter dans l'arborescence fondamentale,
        (iv) le groupe vers lequel on remonte dans l'arborescence fondamentale.

> S:=[
        [p2,0,(u,v,U)->map(t->-t,U),p1],

        [p3,1,(u,v,U)->(map(t->t*exp(2*I*Pi/3),U),map(t->t*exp(-2*I*Pi/3),U)),p1],

        [p4,2,(u,v,U)->map(t->I*t,U),p2],

        [p6,1,(u,v,U)->map(t->t*exp(I*Pi/3),U),p3],

        [cm,0,(u,v,U)->map(t->conjugate(t)*(u+v)/conjugate(u+v),U),p1],

        [pm,0,(u,v,U)->map(t->conjugate(t)*v/conjugate(v),U),p1],

        [pg,0,(u,v,U)->map(t->conjugate(t)*v/conjugate(v)+v/2,U),p1],

        [cmm,2,(u,v,U)->map(t->conjugate(t)*(u+v)/conjugate(u+v),U),p2],

        [pmm,0,(u,v,U)->map(t->conjugate(t)*v/conjugate(v),U),p2],

        [pmg,0,(u,v,U)->map(t->conjugate(t)*v/conjugate(v)+v/2,U),p2],

        [pgg,0,(u,v,U)->map(t->conjugate(t)*v/conjugate(v)+(u+v)/2,U),p2],

        [p31m,1,(u,v,U)->map(t->conjugate(t)*(u+v)/conjugate(u+v),U),p3],

        [p3m1,1,(u,v,U)->map(t->conjugate(t)*v/conjugate(v),U),p3],

        [p4m,2,(u,v,U)->map(t->conjugate(t)*v/conjugate(v),U),p4],

        [p4g,2,(u,v,U)->map(t->conjugate(t)*(v-u)/conjugate(v-u)+(u+v)/2,U),p4],

        [p6m,1,(u,v,U)->map(t->conjugate(t)*v/conjugate(v),U),p6]

]:

On commence par écrire une procédure transfpavages qui
        (i) execute les translations du réseau lorsqu'il s'agit du groupe p1,
        (ii) remonte dans l'arborescence lorsqu'il s'agit d'un autre groupe.

> transfpavages:=proc(type,L,u1,u2,n1,n2)
local M,i,v1,v2;

M:=L;

if type=p1

then return([seq(seq(op(map(U->map(t->t+i*u1+j*u2,U),M)),i=1..n1),j=1..n2)]);

else

        for i from 1 to 16 do

                if type=S[i][1]

                then

                        if S[i][2]=0 then v1:=u1;v2:=u2;

                        elif S[i][2]=1 then v1:=u1;v2:=u1*exp(I*Pi/3);

                        else v1:=u1;v2:=u1*I;

                        end if;

                        return(transfpavages(S[i][4],[op(M),op(map(U->S[i][3](v1,v2,U),M))],v1,v2,n1,n2));

                end if;

        end do;

end if;

end proc:

Enfin, on effectue le tracé de ce que l'on obtient avec transfpavage :

> tracepavage:=proc(type,L,u1,u2,n1,n2)
local M;

M:=transfpavages(type,L,u1,u2,n1,n2);

PLOT(CURVES(op(map(U->map(t->[evalf(Re(t)),evalf(Im(t))],U),M))),AXESSTYLE(NONE),SCALING(CONSTRAINED));

end proc:

> tracepavage(p1,[[0,1.5+I,.5+.8*I,1.5*I,0]],2,2*I+1,3,3);

[Plot]

> tracepavage(p2,[[0,1.5+I,.5+.8*I,1.5*I,0]],2,3*I+1,3,3);

[Plot]

> tracepavage(p3,[[0,1.5+I,.5+.8*I,1.5*I,0]],3,0,3,3);

[Plot]

> tracepavage(p4,[[0,1.5+I,.5+.8*I,1.5*I,0]],3,0,3,3);

[Plot]

> tracepavage(p6,[[0,1.5+I,.5+.8*I,1.5*I,0]],3,0,3,3);

[Plot]

> tracepavage(cm,[[0,1.5+I,.5+.8*I,1.5*I,0]],3,2*I+1,3,3);

[Plot]

> tracepavage(pm,[[0,1.5+I,.5+.8*I,1.5*I,0]],3,1.5*I,3,3);

[Plot]

> tracepavage(pg,[[0,1.5+I,.5+.8*I,1.5*I,0]],3,1.5*I,3,3);

[Plot]

> tracepavage(cmm,[[0,1.5+I,.5+.8*I,1.5*I,0]],3,0,3,3);

[Plot]

> tracepavage(pmm,[[0,1.5+I,.5+.8*I,1.5*I,0]],3.5,3*I-1,3,3);

[Plot]

> tracepavage(pmg,[[0,1.5+I,.5+.8*I,1.5*I,0]],2,6*I+2,6,2);

[Plot]

> tracepavage(pgg,[[0,1.5+I,.5+.8*I,1.5*I,0]],2,6*I+2,6,2);

[Plot]

> tracepavage(p31m,[[0,1.5+I,.5+.8*I,1.5*I,0]],3,0,3,3);

[Plot]

> tracepavage(p3m1,[[0,1.5+I,.5+.8*I,1.5*I,0]],3,0,3,3);

[Plot]

> tracepavage(p4m,[[0,1.5+I,.5+.8*I,1.5*I,0]],3,0,3,3);

[Plot]

> tracepavage(p4g,[[0,1.5+I,.5+.8*I,1.5*I,0]],5,0,3,3);

[Plot]

> tracepavage(p6m,[[0,1.5+I,.5+.8*I,1.5*I,0]],3,0,3,3);

[Plot]