clear all
m=800000; %number of iterations

%Five Transformations
A = 1/2*eye(3);
t(:,1)=[0 ; 0 ; 0];
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];


% Initial point
x(1)=0;
y(1)=0;
z(1)=0;

%Initial point as a vector
v=[0;0;0];

for n=2:m/2
    k=rand;
    % 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);
            x(n)=v(1);
            y(n)=v(2);
            z(n)=v(3);
        end
    end
    hold on
    if mod(n,10000) == 0        %Plotting the pyramid
        u = n/10000;
        plot3(x,y,z,'.','Markersize',4)
        view([0,0])             %Initially we only see one face
        axis([0 1 0 1 0 1]);
        set(gcf, 'color', 'white');
        saveas(gcf, ['Sierpinski_Pyramid',sprintf('%03d',u),'.png'])
        pause(0.000000001)
    end
end

for elev = 0:20                 %Rotate the pyramid to get to the new 3D-view
    u = 41+elev;
    hold on
    plot3(x,y,z,'.','Markersize',4)
    view([elev,floor(elev/2)])
    axis([0 1 0 1 0 1]);
    set(gcf, 'color', 'white');
    saveas(gcf, ['Sierpinski_Pyramid',sprintf('%03d',u),'.png'])
    pause(0.000000001)
end

for n = m/2:m
    k=rand;
    % 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);
            x(n)=v(1);
            y(n)=v(2);
            z(n)=v(3);
        end
    end
    hold on
    if mod(n,10000) == 0 
        u = 62+n/10000;
        plot3(x,y,z,'.','Markersize',4)
        view([20,10])
        axis([0 1 0 1 0 1]);
        set(gcf, 'color', 'white');
        saveas(gcf, ['Sierpinski_Pyramid',sprintf('%03d',u),'.png'])
        pause(0.000000001)
    end
end