 function  envelope(A,t,m)
 % The function 'envelope(A,t,m)' plots the envelope,
 % using 't' rotations and an 'm'-by-'m' partition,
 % and the eigenvalues of a square matrix 'A'.
 %
 % Authors: Panayiotis Psarrakos and Michael Tsatsomeros.

 figure, tic, hold on % The command tic-toc computes the execution time.
 M=moviein(t);
 for j=1:t, cshell(A,m,(2*pi/t)*(j-1)), M(:,j)=getframe; end
 sp(A), hold off, toc, xlabel('Real Axis'), ylabel('Imaginary Axis')

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function cshell(A,N,theta)
 % The function 'cshell(A,N,theta)' plots the shell of a
 % square matrix 'A', with respect to the angle 'theta',
 % using an 'N'-by-'N' partition.

 if nargin<2, k=500; theta=0; end
 [n,mx]=size(A); B=exp(i*theta)*A; R=1.2*norm(A,2);
 xmin=-R; xmax=R; ymin=-R; ymax=R;
 tolx=(xmax-xmin)/N; toly=(ymax-ymin)/N;
 hp=(B+B')/2; shp=(B-B')/2;  % The hermitian part of 'B'.
 [S,T]=eig(hp);  % Solve the eigenproblem of 'hp'.
 % The eigenvalues of the hermitian part 'hp' are sorted
 % in ascending order and the 2 largest eigenvalues 'd1'
 % and 'd2' of 'hp' are found.
 [sa,I]=sort(real(diag(T))); d1=sa(n); d2=sa(n-1);
 % A normalized eigenvector of 'hp' corresponding to 'd1'.
 y=S(:,I(n)); y=y/norm(y,2);
 % Find 'v', that is, the quantity v(B).
 aux=shp*y; v=norm(aux,2)^2;
 % Find 'u', that is, the quantity u(B), but avoid
 % unnecessary computations when u = 0.
 if isreal(B), u=0; else, u=imag(y'*aux); end
 % Compute the values of the two-variable polynomial
 % at the partition points.
 for I=1:N, xx=(xmin+(I-1)*tolx);
     for J=1:N
         yy=(ymin+(J-1)*toly); z=exp(i*theta)*(xx+i*yy);
         xxx=real(z); yyy=imag(z); X(I)=xx; Y(J)=yy;
         Z(I,J)=(d2-xxx)*((d1-xxx)^2+(u-yyy)^2)+(d1-xxx)*(v-u^2);
     end, end
 % Plot the shell (cubic curve).
 contour(X,Y,exp(-i*theta)*Z',[0 0],'b')

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function  sp(A)
 % The function 'sp(A)' plots the eigenvalues of a square
 % matrix 'A' as red +'s.

 [n,m]=size(A); hold on, h=eig(A);
 for j=1:n, plot(real(h(j)),imag(h(j)),'r+'), end
