%%firetower
%
%Original Maple version written by Aage Nielsen.
%Modified MATLAB version by Iver Ottosen.

%%
%The number of cubes in the tower is n.
%The side length of the first cube is l1 and s is the half of that length.
%
%I plot the part of the x and y axis between -c and c, and the 
%part of the z axis between 0 and 2c.

n = 12;

l1 = 1;
s = l1/2;

c = 2.5;

%%
%I define the three matrices on page 49 of the article. 
%The '*' in the notation is omitted.

R = [3/sqrt(10) -1/sqrt(10) 0 0; 1/sqrt(10) 3/sqrt(10) 0 0; 0 0 1 0; 0 0 0 1]

S = [sqrt(10)/4 0 0 0; 0 sqrt(10)/4 0 0; 0 0 sqrt(10)/4 0; 0 0 0 1]

T1 = [1 0 0 0; 0 1 0 0; 0 0 1 l1; 0 0 0 1]

%The 'pause' causes the program to stop and wait for you to press any key.

pause

%%
%Now I calculate the product matrix M.

M = T1*S*R

pause
%%
%I insert the homogeneous coordinates of the eight corners of the first
%cube as rows in a matrix CT.

CT = [s -s 0 1; s s 0 1; -s s 0 1; -s -s 0 1; s -s 1 1; s s 1 1; -s s 1 1; -s -s 1 1];

%The transpose C of the matrix CT have these homogeneous coordinates 
%as column vectors.

C = CT'

pause
%%
%Please ignore the following lines. 
%(Here we read the screen size, open a figure window of appropriate size and 
%give the limits for the coordinate axis.)

scrsz = get(0,'ScreenSize');
sz = scrsz(4)/10;

fig = figure('Position',[sz sz 8*sz 8*sz]);

axis([-c c -c c 0 2*c]);
    
%%
%I now plot the n cubes in the tower. 
%
%D is the matrix obtained from C by deleting the last row. The columns 
%of D is the (ordinary) coordinates for the corners of the cube.
%
%'cube(D)' plots the cube associated with the matrix D.

for i=1:n
    D=C(1:3,1:8);
    cube(D);
    C = M*C;
end

