% This is a program developed to plot the q-numerical range of a matrix
% 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 this program
% if you use it.
%
function y = wqrange(q,A)
%
% Function wqrange(q,A)
%
%
H = (A+A')/2;
G = (A-A')/(2*i);
K = A'*A;
%
% You may increase m to get a better approximation.
%
m = 20;
for r=1:(4*m+1)
              T = cos( (r-1)*pi/(2*m) )*H + sin( (r-1)*pi/(2*m) )*G;
  for s=1:(m+1)
     [U,D] = eig(cos( (s-1)*pi/(2*m) )*T + sin( (s-1)*pi/(2*m) )*K);
     [d,tt] = max(real(diag(D)));
     u = U(:,tt);
     P=real(u'*H*u);
     Q=real(u'*G*u);
     S=real(u'*K*u);
     R=sqrt((1-q^2)*(S-P^2-Q^2));
     P = q*P;
     Q = q*Q;
     jj = (m-1)*r+s;
        for k=1:m
           t = k*pi/m;
           X(k,jj) = P + R*cos(t);
           Y(k,jj) = Q + R*sin(t);
           X(m+k,jj) = P - R*cos(t);
           Y(m+k,jj) = Q - R*sin(t);
           X(2*m+1,jj) = X(1,jj);
           Y(2*m+1,jj) = Y(1,jj);
        end
   end
end
%
%
plot(X,Y);
title('The q-numerical range of A');
xlabel('real axis');
ylabel('imaginary axis');
pause;
grid;
%
% Now, display a sentence and end the program.
%
y = 'You have just seen the q-numerical range of A!';