% function [x1,x2,fail,reig_min] = ric_eig(ham,epp)
%
%  Solves the eigenvalue problem associated with the
%  stabilizing solution (A+R*X stable) of the Riccati equation 
%
%               A'*X + X*A + X*R*X - Q = 0.
%
%  An eigenvalue decomposition is used to obtain the stable
%  invariant subspace of the Hamiltonian matrix, HAM,
%  which contains the above variables in the following format:
%
%       HAM = [A  R; Q  -A'].
%
%  If HAM has no jw axis eigenvalues, there exists n x n matrices
%  X1 and X2 such that [ X1 ; X2 ] spans the n-dimensional
%  stable-invariant subspace of HAM. IF X1 is invertible, then
%  X = X2/X1 satisfies the Riccati equation, and results in A+RX being
%  stable. FAIL is returned with a value of 1 if jw-axis eigenvalues 
%  of HAM are found. An eigenvalue is considered to be purely 
%  imaginary if the magnitude of the real part is less than EPP.
%  The minimum real part of the eigenvalues is returned in REIG_MIN.
%  The default value for EPP is 1e-10.
%
%  This can give incorrect results if HAM is not diagonalizable,
%  in this case it is better to use RIC_SCHR.
%
%  See also: EIG and RIC_SCHR.

function [r1,r2,fail,reig_min] = ric_eig(ham,epp)
  if nargin <= 0
    disp('usage: [x1,x2] = ric_eig(ham,epp)')
    return
  end
  [mtype,mrows,mcols,mnum] = minfo(ham);
  if mtype ~= 'cons'
    error(['RIC_EIG is valid only with CONSTANT matrices'])
    return
  end
  if nargin == 1
    epp = 1e-10;
  end
  np = mrows/2;
  fail = 0;
  [v,d] = eig(ham);
  reig_min = min(abs(real(diag(d))));
%
% check there are no jw-axis eigenvalues in the Hamiltonian matrix
%
  if reig_min < epp
    fail = 1;
  else
    v2 = zeros(2*np,np);
    i2 = 1;
    for i=1:2*np
      if real(d(i,i))<0
        v2(:,i2)=v(:,i);
        i2=i2+1;
      end 
    end
    if i2 == np+1
      r1 = v2(1:np,:);
      r2 = v2((np+1):2*np,:);
    else
      fail = 1;
    end
  end
%
% Copyright MUSYN INC 1991,  All Rights Reserved
