////////////////////////////////////////
// genere et trace la configuration finale d'un tas de sable issu de la configuration initiale vide + 4n^2 grain en l'origine
function sandpile(n,h)
    // initialisation 
    M=h+zeros(2*n+1,2*n+1);
    M(n+1,n+1)=4*n^2;
    // on lance l'algo
    while (max(M)>=4),
        [i,j]=find(M==max(M));
        m=length(i);
        for k=1:1:m
            M(i(k),j(k))=M(i(k),j(k))-4;
            M(i(k)+1,j(k))=M(i(k)+1,j(k))+1; 
            M(i(k),j(k)+1)=M(i(k),j(k)+1)+1; 
            M(i(k)-1,j(k))=M(i(k)-1,j(k))+1; 
            M(i(k),j(k)-1)=M(i(k),j(k)-1)+1;  
        end
    end
    // on trace le résultat
    //f = scf();
    //f.color_map = bonecolormap(1);
    Matplot1(2*M+8,[0 0 1 1])
endfunction

////////////////////////////////////////
// procedure d'ecroulement : entree configuration quelconque, sortie configuration stable
function res=ecroule(M)
    // on lance l'algo
    while (max(M)>=4),
        [i,j]=find(M==max(M));
        m=length(i);
        for k=1:1:m
            M(i(k),j(k))=M(i(k),j(k))-4;
            M(i(k)+1,j(k))=M(i(k)+1,j(k))+1; 
            M(i(k),j(k)+1)=M(i(k),j(k)+1)+1; 
            M(i(k)-1,j(k))=M(i(k)-1,j(k))+1; 
            M(i(k),j(k)-1)=M(i(k),j(k)-1)+1;  
        end
    end
    res=M;
endfunction


function exo2(n,m)
    // condition initiale 
    M=zeros(2*n+1,2*n+1);
    L=sqrt(n);
    C=n+1-L:1:n+1+L;
    M(C,C)=grand(2*L+1,2*L+1,'bin',1,1/2);  
    Matplot1(2*M+8,[0 0 1 1])
    // on itere la procedure suivante m fois
   for k=1:1:m
       aux=find(M==1); // on regarde les sites occupés
       l=length(aux) // il y en a un nombre l
       alea=grand(1,1,'uin',1,l); // on en choisit un au hasard
       indice=aux(alea); //on choisit l'indice correspondant dans aux
       M(indice)=M(indice)+3;
       M=ecroule(M);
       Matplot1(2*M+8,[0 0 1 1])
   end
endfunction

function exo2bis(n,m)
    // condition initiale 
    M=zeros(2*n+1,2*n+1);
    L=sqrt(n);
    C=n+1-L:1:n+1+L;
    M(C,C)=3;  
    Matplot1(2*M+8,[0 0 1 1])
    // on itere la procedure suivante m fois
   for k=1:1:m
       aux=find(M>0); // on regarde les sites occupés
       l=length(aux) // il y en a un nombre l
       alea=grand(1,1,'uin',1,l); // on en choisit un au hasard
       indice=aux(alea); //on choisit l'indice correspondant dans aux
       M(indice)=M(indice)+1;
       M=ecroule(M);
       Matplot1(2*M+8,[0 0 1 1])
   end
endfunction