clear;
 
 N = 25; // taille de la grille (noeuds)
 Naff = 20; // taille affichee
 centre = N/2; // centre de la grille
 L = 2; // largeur de la perturbation (loi empirique)
 
 points_initiaux_x = zeros(N, N); // initialisation des
 points_initiaux_y = zeros(N, N); // coordonnées
 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; // distance au centre
       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; // vecteur unitaire
          direction_y(i,j) = points_initiaux_y(i, j)/d;
       end
    end
 end
 
 // fonction
 
 deff('[e] = profil(r)','e = r*exp(-(r/L)^2)'); // profil empirique
 
 // programme principal
 
 unite = 0; dizaine = 0; // pour le nom de fichier
 for r = 0.1:0.2:centre*1.42 // rayon du front d'onde
    intensite = r^(-1/3); // variation avec la distance
    for i = 1:N
       for j = 1:N
          // e = deplacement du point
          e = intensite*profil(distance(i,j)-r);
          points_finaux_x(i, j) = points_initiaux_x(i, j) + e*direction_x(i, j);
          points_finaux_y(i, j) = points_initiaux_y(i, j) + e*direction_y(i, j);
       end
    end
    // dessin
    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