clear;
 
 N=35; // taille de la grille
 Naff = 30; // 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_f1_x = points_initiaux_x;
 points_f2_x = points_initiaux_x;
 points_f1_y = points_initiaux_y;
 points_f2_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;
       if d == 0 then
          direction_x(i,j) = 0 ; direction_y(i,j) = 0 ;
       else
          direction_x(i,j) = points_initiaux_x(i,j)/d;
          direction_y(i,j) = points_initiaux_y(i,j)/d;
       end
    end
 end
 
 // fonction
 
 deff('[e]=profilP(r)','e=r*exp(-(r/L)^2)'); // profil empirique
 deff('[e]=profilS(r)','e=r*exp(-(r/L)^2)');
 
 // programme principal
 
 unite=0;dizaine=0;centaine=0; // pour le nom de fichier
 for r=0.1:0.2:centre*3.14 // rayon du front d'onde
    intensiteP = r^(-1/3); // variation avec la distance
    intensiteS = %pi/8*r^(-1/2);
    for i=1:N
       for j=1:N
          // e = deplacement du point
          eP = intensiteP*profilP(distance(i,j)-r);
          points_f1_x(i,j) = points_initiaux_x(i,j) + eP*direction_x(i,j);
          points_f1_y(i,j) = points_initiaux_y(i,j) + eP*direction_y(i,j);
          eS = intensiteS*profilS(distance(i,j)-r*.5);
          x = points_f1_x(i,j) ; y = points_f1_y(i,j);
          points_f2_x(i,j) = x*cos(eS) + y*sin(eS); // rotation
          points_f2_y(i,j) = x*sin(eS) - y*cos(eS);
       end
    end
    // dessin
    clf;
    plot2d(0,0,-1,"010"," ",[-Naff/2,-Naff/2,Naff/2,Naff/2]) 
    for i=1:N
       xpoly(points_f2_x(i,:),points_f2_y(i,:));
       xpoly(points_f2_x(:,i),points_f2_y(:,i));
    end
    nom='test\test'+code2str(centaine)+code2str(dizaine)+code2str(unite)+'.gif';
    xs2gif(0,nom,0);
    unite = unite+1;
    if unite == 10 then
       dizaine = dizaine + 1;
       unite = 0 ;
    end
    if dizaine == 10 then
       centaine = centaine + 1;
       dizaine = 0 ;
    end
 end