function menger(n,dx,dy,dz,ox,oy,oz,nb)
% MENGER(M) définit la matrice M représentant l'éponge de Menger
% N est le niveau de récursivité. DX, DY, DZ représentent le décalage en X,
% Y et Z localement. OX, OY et OZ représentent le décalage en X, Y et en Z au 
% niveau N.
global M

if n==0     
    % fin de la récursivité: on met un 1 dans la matrice M pour la
    % case correspondante
    i=dx+ox+1;
    j=dy+oy+1;
    k=dz+oz+1;
    M(i,j,k)=1;
else
    % décalage suuivant le niveau
    ox=ox+nb*dx;  
    oy=oy+nb*dy;
    oz=oz+nb*dz;
    % appel 20 fois récursivement pour les 27 cases moins les 6 cases des
    % milieux des faces et la case au centre (27-6-1=20);
    menger(n-1,0,0,0,ox,oy,oz,nb/3)
    menger(n-1,0,1,0,ox,oy,oz,nb/3)
    menger(n-1,0,2,0,ox,oy,oz,nb/3)
    menger(n-1,1,0,0,ox,oy,oz,nb/3)
    menger(n-1,2,0,0,ox,oy,oz,nb/3)
    menger(n-1,1,2,0,ox,oy,oz,nb/3)
    menger(n-1,2,1,0,ox,oy,oz,nb/3)
    menger(n-1,2,2,0,ox,oy,oz,nb/3)
    menger(n-1,0,0,1,ox,oy,oz,nb/3)
    menger(n-1,0,2,1,ox,oy,oz,nb/3)
    menger(n-1,2,0,1,ox,oy,oz,nb/3)
    menger(n-1,2,2,1,ox,oy,oz,nb/3)
    menger(n-1,0,0,2,ox,oy,oz,nb/3)
    menger(n-1,0,1,2,ox,oy,oz,nb/3)
    menger(n-1,0,2,2,ox,oy,oz,nb/3)
    menger(n-1,1,0,2,ox,oy,oz,nb/3)
    menger(n-1,2,0,2,ox,oy,oz,nb/3)
    menger(n-1,1,2,2,ox,oy,oz,nb/3)
    menger(n-1,2,1,2,ox,oy,oz,nb/3)
    menger(n-1,2,2,2,ox,oy,oz,nb/3)
end
