% This is a program developed to plot the Krein space numerical range 
% 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. 
%
% krein(S,A,m,tol), where S is the hermitian matrix that defines the 
%        indefinite inner product, A is the Krein space operator,
%        4*m^2 is the number of boundary points (x,y,z) on W(H,G,S) that 
%        the program evaluate, and tol > 0 is the lower limit under
%        which the program will compute a point (x/z, y/z) in the 
%        Krein space numerical range.
%
%
function L = krein(K,A,m,d)
%
R = (K*A)/2; 
H = R'+R;
G = (R-R')/i;
%
     [UU,DD] = eig(K);
     [dd,tt] = max(real(diag(DD)));
     uu = UU(:,tt);
     xx = real(uu'*H*uu)/dd;
     yy = real(uu'*G*uu)/dd;
%
% number of iterations and tolerance for the size of the z-coordinate
% are recorded as m and d, respectively
%
for r=1:(4*m+1)
              T = (cos((r-1)*pi/(2*m)))*H + (sin((r-1)*pi/(2*m)))*G;
   for s=1:(2*m)
     [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);
     z1 = real(u'*K*u);  Z(r,s) = z1;
                  if z1 > d,
                      X(r,s)=real(u'*H*u)/z1;
                      Y(r,s)=real(u'*G*u)/z1;
                    else
                      X(r,s) = xx;
                      Y(r,s) = yy;
                    end
   end
end
%
%
      plot(xx,yy,'o',X,Y);
      title('The S-Numerical Range of A');
      grid; 
%