% This is a program for plotting the c-numerical range of A
% for a real n vector c and an nxn complex matrix A by Chi-Kwong Li.
% Please feel free to use or modify it for your research or study.
% It would be nice if you would acknowledge the origin of the program
% if you use it.
%
function y = w(c,a)
%
% Function w(c,A)
%
% Used to plot the convex hull of the c-numerical range of A, which is on
% the inside of the boxes which are drawn. A is a square matrix, and
% size(A) = length(c).
% Uses C.K. Li's algorithm in LMA, 1989, vol 25, pp 285-286.
%
c = -(sort(-c));
H = (a+a')/2;
G = (a-a')/2/sqrt(-1);
%
% Next, the supporting lines of Wc(A) are computed.  A is rotated from
% 0 to pi/2 radians, then the maxima of the edges of Wc(A) are recorded
% to draw them as boxes later.
%
for k = 1:50
   t = k/50*pi/2;
   Ht = cos(t)*H-sin(t)*G;
   Gt = cos(t)*G+sin(t)*H;
   h = real(eig(Ht));
   g = real(eig(Gt));
   a1 = c * sort(h);
   a2 = c * (-sort (-h));
   b1 = c * sort(g);
   b2 = c * (-sort (-g));
   A1 = [a1 a1 a2]';
   B1 = [b1 b2 b2]';
   A2 = [a1 a2 a2]';
   B2 = [b1 b1 b2]';
   x(:,2*k-1) = cos(t)*A1+sin(t)*B1;
   y(:,2*k-1) =-sin(t)*A1+cos(t)*B1;
   x(:,2*k) = cos(t)*A2+sin(t)*B2;
   y(:,2*k) =-sin(t)*A2+cos(t)*B2;
end
%
% Now, the rectangles enclosing Wc(A) are plotted.
%
plot(x,y);
title('The c-Numerical Range of A');
xlabel('real axis');
ylabel('imaginary axis');
grid;
pause;
%
% Now, display a sentence and end the program.
%
y = 'YOU HAVE JUST SEEN THE c-NUMERICAL RANGE OF A';