% This is a program developed to plot the joint numerical range of 3
% hermitian matrices by Chi-Kwong Li.  Please feel free to use or
% modify it for your study or research.  It would be nice if you would
% aknowledge the origin of the program in case you use it.
%
% W(H,G,K,m), where H,G,K should be hermitian matrices of the same
%             size, and the program will evaluate 4m^2  boundary  
%             points of the joint numerical range,
%
% is used to plot the inside convex polytope of the joint umerical range  
% of 3 hermitian matrices H, G, K.
%
function L = wjoint(H,G,K,m)
%
% 
%
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(sin( (s-1)*pi/(2*m) )*T + cos( (s-1)*pi/(2*m) )*K);
     [d1,t1] = max(real(diag(D)));
     u = U(:,t1);
     X(r,s)=real(u'*H*u);
     Y(r,s)=real(u'*G*u);
     Z(r,s)=real(u'*K*u);
     [d2,t2] = max(real(diag(-D)));
     v = U(:,t2);
     X(r,m+1+s)=real(v'*H*v);
     Y(r,m+1+s)=real(v'*G*v);
     Z(r,m+1+s)=real(v'*K*v);
   end
end
%
%
    meshc(X,Y,Z);
    title('The Joint Numerical Range of (H,G,K)');
    xlabel('axis for H');
    ylabel('axis for G');
    zlabel('axis for K');
%