## Copyright (C) 1996, 1998 Auburn University.  All rights reserved.
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## Octave is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
## for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, write to the Free
## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.

## usage: [Xinf,x_ha_err] = hinfsyn_ric(A,BB,C1,d1dot,R,ptol)
##
## forms
##        xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1 BB'];
##        Ha = [A 0*A; -C1'*C1 -A'] - xx;
## and solves associated Riccati equation
## returns error code
##  x_ha_err:
##    0: successful
##    1: Xinf has imaginary eigenvalues
##    2: Hx not Hamiltonian
##    3: Xinf has inf. eigenvalues (numerical overflow)
##    4: Xinf not symmetric
##    5: Xinf not positive definite
##    6: R is singular

function [Xinf, x_ha_err] = hinfsyn_ric (A, BB, C1, d1dot, R, ptol)

  x_ha_err = 0;        # assume success
  Xinf = [];		     # default return value
  n = is_square(A);
  nw = is_square(R);
  if(rank(R) != nw)    x_ha_err = 6;
  else                 # build hamiltonian Ha for X_inf
    xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1, BB'];
    Ha = [A, 0*A; -C1'*C1, -A'] - xx;
    x_ha_err = 0;
    [d, Ha] = balance(Ha);
    [u, s] = schur(Ha, "A");
    rev = real(eig(s));

    if (any(abs(rev) <= ptol))	# eigenvalues near the imaginary axis
      x_ha_err = 1;
    elseif (sum(rev > 0) != sum(rev < 0))
      ## unequal number of positive and negative eigenvalues
      x_ha_err = 2;
    else
      ## compute positive Riccati equation solution
      u = d * u;
      Xinf = u(n+1:2*n,1:n) / u(1:n,1:n);
      if (!all(all(finite(Xinf))))
	x_ha_err = 3;
      elseif (norm(Xinf-Xinf') >= 10*ptol)
	## solution not symmetric
	x_ha_err = 4;
      else
	## positive semidefinite?
	## force symmetry (faster, avoids some convergence problems)
	Xinf = (Xinf + Xinf')/2;
	rev = eig(Xinf);
	if (any(rev <= -ptol))
	  x_ha_err = 5;
	endif
      endif
    endif
  endif
endfunction