% function out = hinfnorm(sys,ttol)
%
%   Calculates the H_infinity norm of stable, SYSTEM
%   matrices, using the Hamiltonian method, or of VARYING
%   matrices, using PKVNORM. The second argument is used
%   for SYSTEM matrices, and is the relative tolerance
%   between the upper and lower bounds for the infinity
%   norm when the search is terminated. There is an optional
%   third argument, for a initial frequency guess, if desired.
%   for SYSTEM matrices, OUT is a 1x3 row vector, with the
%   upper bound, lower bound, and the frequency where the peak
%   occurs. The default value for TTOL is 0.001.
%
%   See also: H2SYN, H2NORM, HINFCHK, HINFSYN, and HINFFI.

function out = hinfnorm(sys,ttol,iiloc)
 if nargin == 0
   disp('usage: out = hinfnorm(sys,ttol)')
   return
 end
 if nargin == 1
   ttol = 0.001;
   iiloc = 5*abs(rand(1,1));
 end
 if nargin == 2
   iiloc = 5*abs(rand(1,1));
 end
 [mtype,mrows,mcols,mnum] = minfo(sys);
 if mtype == 'vary'
   out = pkvnorm(sys);
 elseif mtype == 'cons'
   if nargout == 1
     out = norm(sys);
   else
     disp(['SYSTEM is a constant, norm is ' num2str(norm(sys))])
   end
 else
   [a,b,c,d] = unpck(sys);
   if max(real(eig(a))) >= 0
     error('SYSTEM has closed-right-half plane poles')
     return
   end
   idn = eye(mnum);
   [p,hesa] = hess(a);
   hesb = p' * b;
   hesc = c * p;
   g = hesa \ hesb;
   tfzero = d - hesc*g;
   ntfz = norm(tfzero);
   tfiiloc = (sqrt(-1)*iiloc*idn - hesa) \ hesb;
   nnn = norm(hesc*tfiiloc + d);
   nd = norm(d);
   if nd >= ntfz
     if nd >= nnn
       nrma = nd;
       loc = inf;
     else
       nrma = nnn;
       loc = iiloc;
     end
   else
     if ntfz >= nnn
       nrma = ntfz;
       loc = 0;
     else
       nrma = nnn;
       loc = iiloc;
     end
   end
   if nrma == 0
     nrmt = 1e-8;
   else
     nrmt = (1+ttol)*nrma;
   end
   code = hinfchk(sys,nrmt,1e-9);
   j = 0;

   maxit = 100;
   while j < maxit
     j = j+1;
     if code == -1
       if nargout == 0
         disp(['norm between ' num2str(nrma) ' and ' num2str(nrmt)]);
         disp(['achieved near ' num2str(loc)]);
       else
         out = [nrma nrmt loc];
       end
       j = maxit + 1;
     else
       if length(code) == 1
         probfreq = code;
       else
%        probfreq = code(1:length(code)-1) + 0.5*diff(code);
         probfreq = code(1:length(code)-1) + 0.5*diff(code);
       end
       tmp = 0;
       for jj=1:length(probfreq)
         g = ( (sqrt(-1)*probfreq(jj)) * idn - hesa) \ hesb;
         try = norm(d + hesc*g);
         if try > tmp
           tmp = try;
           loc = probfreq(jj);
         end
       end
       nrma = tmp;
       nrmt = (1+ttol)*nrma;
       code = hinfchk(sys,nrmt,1e-9);
      end
    end
  end
%
% Copyright MUSYN INC 1991,  All Rights Reserved
