function  z = realimagcheck(A, B, eps, all, eigtol)
% Seplambda test: look for a z such that 
%     sigma_min(A - zI) = sigma_min(B - zI) = eps
% If there isn't one, return nan.
% all: 1 if all z's are wanted (may return many duplicates), 0 if one is enough
% eigtol: tolerance used for real, imaginary parts zero (could cause trouble)

if nargin < 4
   all = 0;
end
if nargin < 5
   eigtol = 1e-8;    % very arbitrary
end
nA = length(A); 
nB = length(B);
IA = eye(nA);
IB = eye(nB);
OA = zeros(nA);
OB = zeros(nB);
% we are looking for z = alpha + i*beta for which
%  [      0                 A - (alpha + i*beta)I
%   A' - (alpha - i beta)I          0            ]   has eigenvalue eps, i.e.
%  [    -eps                A - (alpha + i*beta)I
%   A' - (alpha - i beta)I         -eps          ]   has determinant 0, i.e.
%  AA - alpha JJA has eigenvalue i*beta, where
AA = [ A     -eps*IA        % this is Hamiltonian
       eps*IA   -A'  ];   
IIA = eye(2*nA);
IIB = eye(2*nB);
JJA = [ IA   OA
        OA  -IA ];
%  and likewise for B, i.e. BB - alpha JJB has eigenvalue i*beta, where
BB = [ B     -eps*IB        % this is Hamiltonian
       eps*IB   -B'  ];
JJB = [ IB   OB
        OB  -IB ];
%  Thus we are looking for alpha for which AA and BB have the SAME eigenvalue
%  i*beta, i.e. det(K - alpha L) = 0, which is a generalized eigenvalue
%  problem with half the eigenvalues being infinite (because of cancellation
%  in L).  This pencil is known to be regular (not singular) as long as 
%  A and B have no common eigenvalue.  See Section 3 of Gu/Overton.
%  The finite eigenvalues have skew Hamiltonian symmetry around the
%  real axis (i.e., Hamiltonian symmetry if multiplied by i).  Should worry
%  about rounding; ideally would use a special skew-Hamiltonian code.
%
%  Previously we had the following, for which eig(K,L) are exactly the
%  same, but then the notation is not consistent with the paper
%%% K = kron(AA,IIB) - kron(IIA,conj(BB')); 
%%% L = kron(JJA,IIB) - kron(IIA,JJB');
K = kron(IIB,AA) - kron(conj(BB'),IIA);  % transpose, NOT conjugate transpose!!!
L = kron(IIB,JJA) - kron(JJB',IIA);
lambda = eig(K,L);   % can't use "eigs" as K, L are indefinite
% Now we have to check ALL the finite REAL eigenvalues alpha, and for each one, 
% find out if the common eigenvalues of AA - alpha JJ and BB - alpha JJ 
% include an IMAGINARY number i*beta, and furthermore, if so, whether the 
% singular values of A - z I and B - z I for z = alpha + i*beta that
% are, by construction, equal to eps are actually in both cases the MINIMUM 
% singular value - the others are of no interest.

real_eigs = find(abs(imag(lambda)) <= eigtol & abs(real(lambda)) < inf);  
z = nan;       % default if don't find a candidate z
found = 0;
for j=1:length(real_eigs);
   alpha = real(lambda(real_eigs(j)));  % remove imaginary rounding
   AAshift = [ A - alpha*IA     -eps*IA
               eps*IA        -A' + alpha*IA ];   
   BBshift = [ B - alpha*IB     -eps*IB
               eps*IB        -B' + alpha*IB ];   
   lamA = eig(AAshift);
   lamB = eig(BBshift);
   % disp('computed lamA, lamB')
   imag_eigsA = find(abs(real(lamA)) <= eigtol);
   imag_eigsB = find(abs(real(lamB)) <= eigtol);
   if ~isempty(imag_eigsA) & ~isempty(imag_eigsB)
      imag_valsA = imag(lamA(imag_eigsA));   % remove real rounding
      imag_valsB = imag(lamB(imag_eigsB));   % remove real rounding
      diff = imag_valsA*ones(1, length(imag_valsB)) - ... 
       ones(length(imag_valsA), 1)*imag_valsB';
      [indA,indB] = find(abs(diff) < eigtol);
      if ~isempty(indA)  % go through all matching pairs
         for iA = indA'   % index vector must be row vector
            for iB = indB'
               beta = (imag_valsA(iA) + imag_valsB(iB))/2;  % same to rounding
               zcandidate = alpha + beta*i;
               % found a common eigenvalue
               svalsA = svd(A - zcandidate*IA);  % singular value check
               svalsB = svd(B - zcandidate*IB);
               if max([abs(min(svalsA) - eps)  abs(min(svalsB) - eps)]) <= eigtol
                  found = found + 1;
                  z(found) = zcandidate;
                  if ~all
                     return;  % one z is enough
                  end
               else
                  % disp('min singular value not close to eps')
               end
            end
         end
      % disp('after double loop')
      end
   end
end
