%%anilintrans
%
%Original version written by Aage Nielsen
%
%This program makes 3D animations of contractions/expansions 
%and shears and combinations of these. See Lay Section 1.9.
%%
%First I find position and size of the figure window.  
%Please ignore these lines.

scrsz = get(0,'ScreenSize');
sz = scrsz(3)/2; 

%%
%The number of pictures in an animation is n+1.
%I plot the part of the x, y and z axis 
%between c and d.

n = 200;
c = 0;
d = 3;

%%
%Now I make the identity matrix. It is easy to change to 
%another matrix if I or you(?) would like to do so.

K = [1 0 0; 0 1 0; 0 0 1]

%%
%Now I make a matrix for a shear transformation. In fact I 
%make a function of r, which for a given value of r calculates 
%the corresponding shear matrix.

Msh = @(r) [1 r 0; 0 1 0; 0 0 1];

%Now I multiply this matrix on K. Meaning that I make shear 
%transformations of the columns of K.

MSH = @(r) Msh(r)*K;
 
%I show what happens for two values of r. The 'pause' causes
%the program to stop and wait for you to press any key.

pause
MSH(2)
MSH(0)

%I now plot the parallelepiped, spanned by the columns 
%of MSH(0). This is a cube.

pause
fig = figure('Position',[sz 1 sz sz]);

axis([c d c d c d])
parallelepiped(MSH(0));

%%
%Now I make an animation, which shows a series of 
%pictures of parallelepipediums spanned by the columns 
%of MSH(t), where t starts at 0 and ends at 1.

pause
figure(fig)

for s=0:n
    t=s/n;
    clf
    axis([c d c d c d])
    parallelepiped(MSH(t));
    drawnow
end

%%
%Now I make a matrix (actually a matrix function) for a 
%contraction (e<1) or an expansion (e>1).
%In which direction does the contraction/expansion (c/e) work?

Mex = @(e) [1 0 0; 0 e 0; 0 0 1 ];

%I make c/e on the columns of K.

MEX = @(e) Mex(e)*K;

%And show one

pause
MEX(2)
 
%I plot the parallelepiped spanned by the columns of MEX(2)

pause
figure(fig)

clf
axis([c d c d c d]);
parallelepiped(MEX(2));

%And animate

pause
figure(fig)

for s=0:n
    t=2*s/n;
    clf
    axis([c d c d c d])
    parallelepiped(MEX(t));
    drawnow
end

%%
%I now multiply the shear matrix on the c/e matrix

Mse = @(s,e) Msh(s)*Mex(e);

%And multiply it on K

MSE = @(s,e) Mse(s,e)*K;

pause
MSE(1,1)
MSE(0,1)
 
%Now I animate. I can only animate according to one variable 
%(time). Therefore there is an expression in t in both
%variables in MSE. t goes from 1 to 2. For t=1, the shear 
%constant is 0, that is, there is no shear transformation.
%For t=1, the c/e constant is 1, that is, there is no 
%contraction/expansion. All together that means, that the 
%first picture is the cube spanned by the columns of K.
%Find out by what the last picture looks like.

pause
figure(fig)

for s=0:n
    t=s/n+1;
    clf
    axis([c d c d c d])
    parallelepiped(MSE(t-1,t));
    drawnow
end

%%
%Find out - pencil and paper - if the shear and the c/e 
%transformations commute. ( They do not).
%I make the same animation as the one I just made, but 
%changes the order of the two transformations.

%Try to explain from the animations, what it means, that 
%the shear and the c/e transformations do not commute.

Mes = @(s,e) Mex(e)*Msh(s);
MES = @(s,e) Mes(s,e)*K;

pause
MES(1,1)
MES(0,2)

pause
figure(fig)

for s=0:n
    t=s/n+1;
    clf
    axis([c d c d c d])
    parallelepiped(MES(t-1,t));
    drawnow
end
