Source
 clear;
 
 N=25; // taille de la grille
 Naff = 20; // taille affichee
 centre = N/2; // centre de la grille
 L=2; // largeur (loi empirique)
 
 points_initiaux_x = zeros(N,N); // initialisation
 points_initiaux_y = zeros(N,N);
 for i=1:N
    for j=1:N
       points_initiaux_x(i,j) = i-centre;
       points_initiaux_y(i,j) = j-centre;
    end
 end
 points_finaux_x = points_initiaux_x;
 points_finaux_y = points_initiaux_y;
 for i=1:N
    for j=1:N
       d = sqrt(points_initiaux_x(i,j)^2 + points_initiaux_y(i,j)^2);
       distance(i,j) = d;
    end
 end
 
 // fonctions
 
 deff('[e]=profil(r)','e=r*exp(-(r/L)^2)');
 
 // programme principal
 
 unite=0;dizaine=0;
 for r=0.1:0.2:centre*1.42 // rayon du front d'onde
    intensite = %pi/8*r^(-1/2); // variation avec la distance
    for i=1:N
       for j=1:N
          e = intensite*profil(distance(i,j)-r);
          x = points_initiaux_x(i,j) ; y = points_initiaux_y(i,j);
          points_finaux_x(i,j) = x*cos(e) + y*sin(e); // rotation
          points_finaux_y(i,j) = x*sin(e) - y*cos(e);
       end
    end
    clf;
    plot2d(0,0,-1,"010"," ",[-Naff/2,-Naff/2,Naff/2,Naff/2]) 
    for i=1:N
       xpoly(points_finaux_x(i,:),points_finaux_y(i,:));
       xpoly(points_finaux_x(:,i),points_finaux_y(:,i));
    end
    nom='test\test'+code2str(dizaine)+code2str(unite)+'.gif';
    xs2gif(0,nom,1);
    unite = unite+1;
    if unite == 10 then
       dizaine = dizaine + 1;
       unite = 0 ;
    end
 end