% JMS - 5 juillet 2003 - Spirale de Fraser 
% Voir : Les illusions des sens - Pour la Science avril/juin 2003
%        Page 38 (illusions géométriques)
% ===============================================================

% unité
u = 15;      

% rayon - la valeur 1.17 est une approximation de la solution
% de l'équation : x - 1 = Pi/20 * sqrt(x) qui permet d'obtenir des 
% dalles proches d'un carré.
def r(expr x) =
    mexp(x * mlog(1.17))
enddef;

% épaisseur des couronnes blanc/noir
def e(expr x) =
    1 + 0.4 * r(x)
enddef;

% point du dallage
def p(expr x,y) =
    u*r(y)*(cosd(x*9),sind(x*9))
enddef;

% point inférieur d'une couronne
def ti(expr x,y) =
    (u*r(y) - (e(y)/2))*(cosd(x*9),sind(x*9))
enddef;
% point inférieur intermédiaire d'une couronne
def tii(expr x,y) =
    (u*r(y) - (e(y)/4))*(cosd(x*9),sind(x*9))
enddef;
% point supérieur d'une couronne
def ts(expr x,y) =
    (u*r(y) + (e(y)/2))*(cosd(x*9),sind(x*9))
enddef;
% point supérieur intermédiaire d'une couronne
def tsi(expr x,y) =
    (u*r(y) + (e(y)/4))*(cosd(x*9),sind(x*9))
enddef;

% cadre de la figure
path cadre;
cadre = ((-1,-1)--(-1,1)--(1,1)--(1,-1)--cycle) scaled (u*r(18));


beginfig(1);

 fill cadre withcolor red+green;

 for i:=9 upto 21:
    for j:=0 upto 40:
     fill p(j,i+0.05) -- p(j+0.45,i+0.5) -- 
          p(j,i+0.95) -- p(j-0.45,i+0.5) -- cycle
          withcolor red;
    endfor;
 endfor;

 for i:=9 upto 21:
    for j:=0 step 2 until 40:
     jj := if i mod 2 = 0: j else: j-1 fi;
     fill p(jj+0.5,i-0.45) -- p(jj+0.95,i) -- 
          p(jj+0.5,i+0.45) -- p(jj+0.05,i) -- cycle
          withcolor .6 * blue;
    endfor;
    draw fullcircle scaled (2*u*r(i))
         withpen pencircle scaled (e(i))
         withcolor white;
    for j:=0 step 2 until 40:
     jj := if i mod 2 = 0: j+1 else: j fi;
     fill (ti(jj-0.7,i) .. ti(jj-0.35,i) .. ti(jj,i)) --
          (ti(jj,i) .. tii(jj+0.5,i) .. ts(jj+1.7,i)) --
          (ts(jj+1.7,i) .. ts(jj+1.35,i) .. ts(jj+1,i)) --
          (ts(jj+1,i) .. tsi(jj+0.5,i) .. ti(jj-0.7,i)) -- cycle
         withcolor black;
    endfor;
 endfor;

 clip currentpicture to cadre;

endfig;

beginfig(2);
 
 fill cadre withcolor red+green;

 for i:=9 upto 21:
    for j:=0 upto 40:
     fill p(j,i+0.05) -- p(j+0.45,i+0.5) -- 
          p(j,i+0.95) -- p(j-0.45,i+0.5) -- cycle
          withcolor red;
    endfor;
 endfor;

 for i:=9 upto 21:
    for j:=0 step 2 until 40:
     jj := if i mod 2 = 0: j else: j-1 fi;
     fill p(jj+0.5,i-0.45) -- p(jj+0.95,i) -- 
          p(jj+0.5,i+0.45) -- p(jj+0.05,i) -- cycle
          withcolor .6 * blue;
    endfor;
    draw fullcircle scaled (2*u*r(i))
         withpen pencircle scaled (e(i))
         withcolor white;
    for j:=0 step 2 until 40:
     jj := if i mod 2 = 0: j else: j-1 fi;
     fill (ti(jj-0.7,i) .. ti(jj-0.35,i) .. ti(jj,i)) --
          (ti(jj,i) .. tii(jj+0.5,i) .. ts(jj+1.7,i)) --
          (ts(jj+1.7,i) .. ts(jj+1.35,i) .. ts(jj+1,i)) --
          (ts(jj+1,i) .. tsi(jj+0.5,i) .. ti(jj-0.7,i)) -- cycle
         withcolor black;
    endfor;
 endfor;
 
 clip currentpicture to cadre;
 
endfig;

end