% Matlab code for a beautiful Sierpinski pyramid
% just type "pyramid" in the Matlab command window to run it
clear all
iterations=40000; %the number of iterations

% recall, the transformations are of the form A*[x,y]+t
% and, there are five such transformations ...
A = 1/2*eye(3);         % How fortunate! they all have the same A.
t(:,1)=[0 ; 0 ; 0];     % Five different ts.
t(:,2)=[0 ; 1/2 ; 0];
t(:,3)=[1/2 ; 0; 0];
t(:,4)=[1/2 ; 1/2; 0];
t(:,5) = [1/4 ; 1/4; 1/2];


% 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 five transformations to get a new point
  for i = 1:5
      if (k > (i-1)/5 && k < i/5)  % Each transformation has equal probability
          v = 1/2*v+t(:,i);  % Do the transformation
          % now, store this new point
          x(n)=v(1);
          y(n)=v(2);
          z(n)=v(3);
      end
     
  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');