% Matlab code for a beautiful Menger sponge
% just type "sponge" in the Matlab command window to run it
clear all
iterations=50000; %the number of iterations

% recall, the transformations are of the form A*[x,y]+t
% and, there are twenty such transformations ...
A = 1/3*eye(3);         % How fortunate! they all have the same A.
t(:,1)=[0 ; 0 ; 0];     % Alas, the ts are all different.
t(:,2)=[0 ; 0 ; 1/3];
t(:,3)=[0 ; 0; 2/3];
t(:,4)=[0 ; 1/3; 0];
t(:,5) = [0 ; 1/3; 2/3];
t(:,6) = [0 ; 2/3; 0];
t(:,7) = [0 ; 2/3; 1/3];
t(:,8) = [0 ; 2/3; 2/3];
t(:,9) = [1/3 ; 0; 0];
t(:,10) = [1/3 ; 0; 2/3];
t(:,11) = [1/3 ; 2/3; 0];
t(:,12) = [1/3 ; 2/3; 2/3];
t(:,13) = [2/3 ; 0; 0];
t(:,14) = [2/3 ; 0; 1/3];
t(:,15) = [2/3 ; 0; 2/3];
t(:,16) = [2/3 ; 1/3; 0];
t(:,17) = [2/3 ; 1/3; 2/3];
t(:,18) = [2/3 ; 2/3; 0];
t(:,19) = [2/3 ; 2/3; 1/3];
t(:,20) = [2/3 ; 2/3; 2/3];

% the initial point
x(1)=0;
y(1)=0;
z(1)=0;

% but, let's write the (x,y,z) points as a vector, v
v=[0;0;0];

% Now begin the iterations
for n=2:iterations
  % choose a random number, k, between 0 and 1
  k=rand;
  % depending on your random number ...
  % do one of the twenty transformations to get a new point
  for i = 1:20
      if k > (i-1)/20   % Each transformation has equal probability
         if k < i/20
             v = 1/3*v+t(:,i);  % Do the transformation
         end
      end
      % now, store this new point
      x(n)=v(1);
      y(n)=v(2);
      z(n)=v(3);
  end
end

%now, let's plot all those (x,y,z) points that we just computed!
opengl software %this is a fix to ensure that Matlab won't crash :)
hold off
plot3(x,y,z,'.','Markersize',4)
view([0,0])
axis('equal');
box on;
set(gcf, 'color', 'white');