% Matlab code RayleighContour.m
% For "Applied Numerical Linear Algebra",  Figure 5.1
% Written by James Demmel, Jun  2, 1997
%
% Plots contours of Rayleigh quotient rho(u,A) on the unit sphere
%
% Input
%    A = any 3-by-3 symmetric matrix
%
% Output
%    Plot of contours of rho(u,A)
%
[u,d]=eig(A);
ds=-sort(-diag(d));
if ds(1)==ds(3),
% A is a multiple of the identity
  disp(['rho(u,A) = ',num2str(ds(1)),' is constant.'])
else 
% choose contour values CON to plot
% normalize spectrum to [1 lam 0]
  ds=ds-ds(3);
  ds=ds/ds(1);
  lam = ds(2);
  if (lam < 1e-13) 
     lam = 1e-13;
  end
  if lam == 1,
     CON = [0 .001 .01 (.1:.1:.9) .99 .999 1];
  else 
     omlam=1-lam;
     CON = [ 0 lam/100 (lam/10:lam/5:lam*9/10) lam ...
             (lam+omlam/10:omlam/5:lam+omlam*9/10) 1-omlam/100 1 ];
  end
  theta=(0:.025:2*pi);
  C = cos(theta);
  S = sin(theta);
  hold off, clg
  for subplt=1:2
  subplot(1,2,subplt)
  for c = CON,
      if (c <= lam)
         x=sqrt(c)*S;
         y=(sqrt(c)/sqrt(lam))*C;
      else
         t=acos(sqrt(lam/c));
         dt = (pi-2*t)/1000;
         TT = [(t:dt:pi-t),(pi+t:dt:2*pi-t)];
         SS = sin(TT);
         CC = cos(TT);
         x=sqrt(c)*SS;
         y=(sqrt(c)/sqrt(lam))*CC;
      end
      z=sqrt(1-x.^2-y.^2);
      R=find(imag(z)==0);
      x=x(R);y=y(R);z=z(R);
      if (length(x)>0),
         rotateP = u*[x;y;z];
         xP=rotateP(1,:);
         yP=rotateP(2,:);
         zP=rotateP(3,:);
         rotateM = u*[x;y;-z];
         xM=rotateM(1,:);
         yM=rotateM(2,:);
         zM=rotateM(3,:);
         x = [xP,xM];
         y = [yP,yM];
         z = [zP,zM];
         if length(x) > 0 & subplt == 2,
            vue = view;
            vue = vue(3,1:3);
            HS = vue*[x;y;z];
            front = find(HS >= 0);
            x = x(front);
            y = y(front);
            z = z(front);
         end
         if (c < lam)
%           Color values less than lam green
            plot3(x,y,z,'g.')
            axis('square')
            hold on
         elseif (c == lam)
%           Color values less than lam white
            plot3(x,y,z,'.w')
            axis('square')
            hold on
         else
%           Color values greater than lam red
            plot3(x,y,z,'r.')
            axis('square')
            hold on
         end
      end
   end
   xlabel('x = e_1')
   ylabel('y = e_2')
   zlabel('z = e_3')
%  title('Contour of Rayleigh Quotient. Red = max, Green = min')
   if subplt == 1,
      grid
      title('Transparent')
   else
      title('Opaque')
   end
   axis([-1 1 -1 1 -1 1])
   end
end