function BBplot = wberell2(A,k,gr,choice,color)

% This function plots the field of values (FOV) boundary from k
% eigenanalyses using boundary point compressions for A to fill in the 
% boundary gaps with ellipses
% Input :  A   a real or complex n by n matrix
%          k   an integer >= 2 (number of eigenanalyses performed)
%          gr  if gr = 2, we plot all internal 2 by 2 compression FOV boundary
%              ellipses and all computed FOV boundary points (in red)
%          choice = 1, 2, 3, or 4 :  (default = 1)
%             allows to choose between 4 strategies for new rotation angles
% Output:  plot of boundary curve of FOV(A) and BBplot, the approx boundry pt 
%          coordinates in RR^2
% Experimental version, (c) Frank Uhlig, 2012

if nargin == 1, k = 14; gr = 0; choice = 1; color = 'blue';end
if nargin == 2, gr = 0; choice = 1; color = 'blue'; end
if nargin == 3, choice = 1; color = 'blue'; end
if nargin == 4, color = 'blue'; end
kryln = 200; n = length(A); BB = [];
if norm(A*A' - A'*A,1) < norm(A,1) * 10^-12   % extra care to deal with normal matrices:
  'A is  normal', BB = eig(A).';  BB = [real(BB);imag(BB)];
  try
    KK = convhull(BB(1,:),BB(2,:)); BBplot = BB(:,KK); % compute convex hull and plot FOV
    plot(BBplot(1,:),BBplot(2,:),'Color',color); axis equal, grid on, 
  catch
    plot(BB(1,:),BB(2,:),'Color',color); axis equal, grid on, BBplot = BB; end 
  title('Field of values boundaries via 2 by 2 compressions (blue) and eigenanalyses-only (red)','FontSize',15); hold off,
  return, end                               % and quit

H = (A+A')/2; K = (A-A')/(2i); 
[Rj Tj] = eigenpts(H,n,kryln);       % max real axis extension FOV pts
R(:,1) = Rj; T(:,1) = Tj;
[Rj Tj] = eigenpts(K,n,kryln);       % max imaginary axis extension FOV pts
R(:,2) = Rj; T(:,2) = Tj;
X = [diag([R,T]'*A*[R,T]).';0, pi/2, pi, 3*pi/2; R, T]; 
X = [X,X(:,1)]; X(2,end) = 2*pi;     % first 4 FOV edge pts via 2 eigenanalyses

for j = 3:k          % refine FOV boundary points
  s2X = size(X,2); lenEj = [];    
% Now : choose either lenEj block below:
  if choice == 1     % default ( = fastest) choice
    for jj = 1:s2X-1 % reduce max distance of tangent intersect from FOV secants
      lenEj = [lenEj, norm(X(2,jj) - X(2,jj+1))*norm(X(1,jj) - X(1,jj+1))]; 
      if jj == 1, titel = 'Use products of angle and distance differences'; 
    end, end 
  elseif choice == 2 % choose to reduce largest tangents/secant triangle area
    for jj = 1:s2X-1  
      njj = [cos(X(2,jj)), sin(X(2,jj))];       % tangent normals as rows
      njjp1 = [cos(X(2,jj+1)), sin(X(2,jj+1))];
      Thnjj = [njj;njjp1];                      % tangent normals matrix
      xjj = [real(X(1,jj));imag(X(1,jj))];      % field of values boundary pts
      xjjp1 = [real(X(1,jj+1));imag(X(1,jj+1))]; 
      xy = Thnjj\diag(Thnjj * [xjj,xjjp1]);     % tangent intersection pt
      xjjxjjp1 = xjjp1 - xjj;                   % secant vector between edge pts
      sjxy = xy - xjj;                          % tangent intersection to FOV pt
      Atjjp1 = abs(det([xjjxjjp1,sjxy])/2);     % tangent triangle area
      lenEj = [lenEj,Atjjp1]; end, 
    titel = 'Use areas of tangent/secant triangles'; 
  elseif choice == 3 % choose to reduce largest angular spread
    for jj = 1:s2X-1                 % determine angular spreads
       lenEj = [lenEj, norm(X(2,jj) - X(2,jj+1))]; end
    titel = 'Use angular spreads';
  elseif choice == 4 % choose to reduce largest distance spread
    for jj = 1:s2X-1                 % determine distance spreads
       lenEj = [lenEj, norm(X(1,jj) - X(1,jj+1))]; end
    titel = 'Use distance spreads';
  end

  [~, in] = max(lenEj); cj = in(1); thj = (X(2,cj) + X(2,cj+1))/2; % new theta
  if thj < pi, thjp = thj + pi; else thjp = thj - pi; end        % antipode
  Ath = cos(thj)*H + sin(thj)*K;     % new test matrix
  [Rcj Tcj] = eigenpts(Ath,n,kryln); % two new FOV boundary pts
  X = [X(:,1:cj),[Rcj'*A*Rcj;thj;Rcj],X(:,cj+1:end)]; % fit Rcj into X data
  ip = find(X(2,:) > thjp); ip = ip(1)-1;             % fit Tcj into X
  X = [X(:,1:ip),[Tcj'*A*Tcj;thjp;Tcj],X(:,ip+1:end)]; end

for j = 1:size(X,2)-1                % for jth subarc ellipse fill-in:
  xx = X(3:end,j); yy = X(3:end,j+1); 
  v = (yy-dot(xx,yy)*xx); 
  if norm(v) > 10^-14                % to weed out repeated corners
    v = v/norm(v);
    T = [xx,v]; Axv = T'*A*T;        % 2 by 2 compression of for xx and v
    [center major minor alpha] = matrixFOV(Axv); 
    B = plotEllipse2pts(X(1,j),X(1,j+1),center, major, minor, alpha,'black',gr); 
    BB = [BB,B]; hold on, end, end
if gr == 2,
  plot(real(X(1,:)),imag(X(1,:)),'sr'), hold on, end % plot all computed FOV bdry pts

k1 = isfinite(BB(1,:)) == 0; BB(:,k1) = [];  % delete non-numberic data in BB  
k2 = isfinite(BB(2,:)) == 0; BB(:,k2) = [];  
KK = convhull(BB(1,:),BB(2,:)); BBplot = BB(:,KK); % compute convex hull and plot FOV
plot(BBplot(1,:),BBplot(2,:),'Color',color); axis equal, grid on, 
title('Field of values boundaries via 2 by 2 compressions and eigenanalyses','FontSize',15); hold off,
hold off, 

function [center, majorAxis, minorAxis, alpha] = matrixFOV(A)
% This function takes a 2x2 matrix A (often a compression), and finds the
% parameters for the ellipse that represents its field of values. 
% Input : A         - a 2x2 matrix
% Output: alpha     - angle of rotation wrt. real axis in CC
%         majorAxis - half of length of major axis (real value)
%         minorAxis - half of length of minor axis
%         center    - center of ellipse            (given as a complex number)
 center = trace(A)/2; A0 = A - (trace(A)/2)*eye(2); % center FOV(A) at zero
 alpha = real(log(sqrt(-det(A0)/abs(det(A0))))/1i); % find angle of main axis
 majorAxis = sqrt(norm(A0,'fro')^2 + 2*abs(det(A0)))/2; % major axis length
 minorAxis = sqrt(norm(A0,'fro')^2 - 2*abs(det(A0)))/2; % minor axis length
 
  function B = plotEllipse2pts(xj,xjp1,center,a,b,alpha,color,gr)
% This function plots an ellipse given its center, angle of rotation, and 
% major and minor axis lengths.
% Input :  xj        - jth point on FOV boundary
%          xjp1      - j+1st point on FOV boundary
%          center    - center of ellipse  (given as a complex number)
%          a		 - half of length of major axis  (real value)
%          b		 - half of length of minor axis  (real value)
%          alpha     - angle of rotation  (counter-clockwise)        
%          color     - color to plot ellipse; (use 'blue' for example)
%          gr        - if gr is equal to 2: plot full compression ellipses
% Output:  B real coordinates of ellipse pts between xj and xjp1      
ric = [real(center);imag(center)]; % convert to R^2 data :
xj = [real(xj);imag(xj)]; xjp1 = [real(xjp1);imag(xjp1)];   
n1 = [cos(alpha),sin(alpha)];      % use trailing direction from xjp1 towards xj for n1
if n1*(xj - xjp1) < 0, n1 = -n1; end, n2 = [-n1(2),n1(1)]; % major axis and normal vectors
A22 = [n1;n2]; b2j = [n1*xj;n2*ric]; b2jp1 = [n1*xjp1;n2*ric]; 
xlj = A22\b2j; xljp1 = A22\b2jp1;  % projections of FOV pts on major ellipse axis
xcj = xlj + (a^2-norm(xlj-ric)^2)^.5 / norm(xj-xlj) * (xj-xlj); % circular angle pts
xcjp1 = xljp1 + (a^2 - norm(xljp1 - ric)^2)^.5 / norm(xjp1 - xljp1) * (xjp1-xljp1);
ths = acos(n1*(xcj - ric) / (norm(xcj-ric)*norm(n1)));      % sustaining circle angles 
the = acos(n1*(xcjp1 - ric) / (norm(xcjp1-ric)*norm(n1))); 
thuse = sort([ths,the]); ths = thuse(1); the = thuse(2);    % use angles sorted by size
theta = linspace(ths,the,ceil(40*norm(xj - xjp1)/(a+b)+4)); % range of partial ellipse plot  
X(1,:) = a*cos(theta); X(2,:) = b*sin(theta);          % Parametric equations data 
A=real(A22'*X); B = [A(1,:)+ric(1);A(2,:)+ric(2)];     % rotated and shifted pt data
if gr == 2,                        % Plot all compression ellipses if desired
  thall = 0:0.1:2*pi;
  Xall(1,:) = a*cos(thall); Xall(2,:) = b*sin(thall);  % Parametric equations data 
  Aall=real(A22'*Xall); Ball = [Aall(1,:)+ric(1);Aall(2,:)+ric(2)]; % rotated and shifted pts  
  plot(Ball(1,:),Ball(2,:),'Color', color); hold on
  axis equal, grid, hold on, end

function [r t] = eigenpts(A,n,kryln)
% Finds the extreme FOV point generating complex eigenvectors for A
  if n < kryln                    % for small n < kryln use QR algorithm
    [V,D] = eig(A); 
    r = V(:,n); t = V(:,1);       % (use max and min evalue/vector)
  else                            % for large n use Krylov methods
    opts.disp = 0; opts.maxit = 1600;
    [r D] = eigs(A,1,'Lr',opts);  % use eigs
    [t,D] = eigs(A,1,'Sr',opts);
  end
