function  shell(a,k,xx)
% Plots the shell (and spectrum) of an 
% almost skew-symmetric matrix `a'
% `k' is the number of plotted points
% `xx' is the character string passed on to spectrum plotting
% Author: Panos Psarrakos and Michael Tsatsomeros

if nargin<2, k=300; xx='b+'; end
[n,mx]=size(a);
sp=(a+a')/2; sk=a-sp; % symmetric and skew parts of a

if rank(sp)==1
  % (r,y) is the dominant eigenpair of sp; y is automatically normalized by matlab!
  [S,T]=eig(sp); [r,I]=max(max(abs(T)));  
  y=S(:,I); sr=sign(T(I,I));
  'The variance of this almost skew-symmetric matrix is'
  v=norm(sk*y,2)^2 % the variance
  'and the nonzero eigenvalue of its symmetric part is'
  d=max(eig(sp)) % the delta
      h=r/k;  x=0;
         for t=1:k,                   
            x=x+h; 
            q=(r-x)*((v/x)+x-r);    
            rt=roots([1 0 -q]);
            if q>0
              aa=x+i*rt(1);  
              bb=x+i*rt(2);  
              plot(sr*real(aa),imag(aa),'r.')
              hold on
              plot(sr*real(bb),imag(bb),'r.')
              hold on      
            end   
            xlabel('Real   Axis')  
            ylabel('Imaginary    Axis')  
         end
         axis([-1 r+1 -n n])
         spectrum(a,xx)
else
' Matrix is not almost skew-symmetric'
end