function menger_sponge(num_iteration)
    close all;
    clf;
    
    % set number of iterations
    n = num_iteration;
    
    % setting up matrix for the plot
    N = 3^n;
    M = ones(N,N,N);
    
    % setting up empty vertice and faces vectors
    global v f
    v = [];
    f = [];
    
    % recursively remove cubes to form the menger sponge
    M = remove_cubes(M);
    
    % generate faces and vertices for graphing
    plot3d(M,N);
    
    % graphing
    patch('Faces', f, 'Vertices', v, 'FaceColor', 'g', 'FaceLighting', 'flat', 'EdgeColor', 'none');
    axis([0 N 0 N 0 N]);
    axis equal;
    axis off;
    light('Position',[N/2 N/2 N/2],'Style','local');
    light('Position',[N N N],'Style','local');
    
    % animation
%     axis vis3d
%     for i=1:180;
%        saveas(gcf,['menger', sprintf('%i',n), sprintf('%03d',i), '.png'])
%        camorbit(2,0,'data',[0 0 1])
%        drawnow
%     end
end

function new_m = remove_cubes(m)
% if m is 1x1x1, then stop recursion
% else, separate m in a 3x3x3 fashion and remove the middle 7 cubes
% and then recurse on the 20 cubes that are left
    s = size(m,1);
    if s == 1
        new_m = m;
    else
        new_m = m;
        len = s/3; % length of one 'cube'
        % a,b,c are meta-coordinates for the 27 'cubes'
        a = 1:len;
        b = (len+1):(2*len);
        c = (2*len+1):(3*len);
        
        A = [a;b;c]';
        for i_c = A
            for j_c = A
                for k_c = A
                    i = i_c';
                    j = j_c';
                    k = k_c';
                    if needs_to_be_removed(i,j,k,a,b,c)
                        new_m(i,j,k) = 0;
                    else
                        new_m(i,j,k) = remove_cubes(m(i,j,k));
                    end
                end
            end
        end             
    end
end

function bool = needs_to_be_removed(i,j,k,a,b,c)
% determines whether a cube needs to be removed
% based on its coordinates, returns true for the middle 7 cubes
    f1 = isequal([i j k],[b b a]) || isequal([i j k],[a b b]);
    f2 = f1 || isequal([i j k],[b a b]) || isequal([i j k],[b b c]);
    f3 = f2 || isequal([i j k],[c b b]) || isequal([i j k],[b c b]);
    bool = f3 || isequal([i j k],[b b b]);
end

function plot3d(M,N)
% plots the matrix M in voxel form
    for x=1:N
        for y=1:N
            for z=1:N
                if M(x,y,z) == 1
                    plotcube(x,y,z);
                end
            end
        end
    end
end

function plotcube(x,y,z)
% converts a coordinate into a voxel and plots it
    vertices = [x y z;
                x y z-1;
                x y-1 z;
                x-1 y z;
                x y-1 z-1;
                x-1 y-1 z;
                x-1 y z-1;
                x-1 y-1 z-1];
    faces = [1 2 7 4;
             1 3 5 2;
             1 4 6 3;
             6 3 5 8;
             6 8 7 4;
             5 8 7 2];
    global v f
    f = [f; faces+size(v, 1)];
    v = [v; vertices];
end