function nmhinf = normhinf(Z1,Z2,Z3,Z4,Z5)
% [NMHINF] = NORMHINF(A,B,C,D,tol) or
% [NMHINF] = NORMHINF(SS_,tol) computes the HINF norm of the given 
%    state-space realization. Implemented here is a binary search 
%    algorithm of imaginary axis eigenvalue(s) of the Hamiltonian
%
%                        -1           -1
%       H(gam) = | A + BR  D'C     -BR  B'        |
%                |       -1                -1     |
%                |C'(I+DR  D')C    -(A + BR  D'C)'|
%
%    where R = gam^2 I - D'D > 0. HINF norm equals to "gam" when H 
%    has imaginary axis eigenvalue(s).
%
%    Inital guesses of the HINF norm upper/lower bounds are
%
%       Upper Bound: max_sigma(D) + 2*sum(Hankel SV(G))
%       Lower Bound: max{max_sigma(D), max_Hankel SV(G)}.
%
%    The search algorithm stops when two adjacent "gam's" have relative
%    error less than "tol". If no "tol" provided, tol = 0.001.

% R. Y. Chiang & M. G. Safonov 8/91
% Copyright (c) 1991 by the MathWorks, Inc.
% All Rights Reserved.

nag1 = nargin;

if nag1 == 5 | nag1 == 2
   inargs='(a,b,c,d,tol)';
   eval(mkargs(inargs,nargin,'ss'))
end

if nag1 == 4 | nag1 == 1
  inargs='(a,b,c,d)';
  eval(mkargs(inargs,nargin,'ss'))
  tol = 0.001;
end

[rd,cd] = size(d);
dtd = d'*d;
hsv = hksv(a,b,c);

n_l = max([max(svd(d)),max(hsv)]);
n_h = max(svd(d)) + 2*sum(hsv);

while (n_h-n_l) > 2*tol*n_l
    gam = (n_l+n_h)/2;
    r = gam*gam*eye(cd) - dtd;
    ir = inv(r);
    Ham = [a+b*ir*d'*c  -b*ir*b';
           c'*(eye(rd)+d*ir*d')*c -(a+b*ir*d'*c)'];
    lam = abs(real(eig(Ham)));
    ind = find(lam < sqrt(eps));
    [rind,cind] = size(ind);
    if rind == 0     % no imag. lam
       n_h = gam;
    else             % exists imag. lam
       n_l = gam;
    end
end

nmhinf = gam;

%
% ----- End of NORMHINF.M  % RYC/MGS %


