function [gpeak,fpeak] = norminf(D,tol,fBand,NeedStable)
% Compute the peak gain GPEAK of the frequency response
%                                   -1
%             G (s) = D + C (sE - A)  B .
%
%   The norm is finite if and only if (A,E) has no eigenvalue on the
%   imaginary axis.  TOL is the desired relative accuracy on GPEAK,
%   and FPEAK is the frequency such that:
%
%           || G ( j * PEAKF ) ||  =  GPEAK     (continuous)
%
%                   j*FPEAK*Ts
%           || G ( e           ) ||  =  GPEAK   (discrete)

%  Author(s):  P. Gahinet
%   Copyright 1986-2011 The MathWorks, Inc.

%    Based on the algorithm described in
%        Bruisma, N.A., and M. Steinbuch, A Fast Algorithm to Compute
%        the Hinfinity-Norm of a Transfer Function Matrix, Syst. Contr.
%        Letters, 14 (1990), pp. 287-293.
Ts = abs(D.Ts);
ni = nargin;
if ni<3
   fBand = [];
end
if ni<4
   NeedStable = false;
end

% Determine if internal delays are equivalent to input+output delays,
% in which case they can be ignored
hasInternalDelay = norm_hasInternalDelay(D);
if hasInternalDelay
   if Ts==0
      throw(ltipack.utNoDelaySupport('norm',0,'internal'))
   else
      D = elimDelay(D,[],[],D.Delay.Internal);
   end
end

% Check stability before simplifying and scaling to avoid inconsistencies
if NeedStable && ~isstable(D)   
   gpeak = Inf;  fpeak = NaN;  return
end

% Simplify and scale realization
D = sminreal(D);
[a,b,c,d,~,e] = getABCDE(D);
if ~D.Scaled
   % Scale for freq. response accuracy
   [a,b,c,e] = ltipack.xscale(a,b,c,d,e,Ts,'Warn',false);
end

% In CT, eliminate cancellations at Inf because frequency response can't be
% reliably evaluated near s=Inf when E is singular
if Ts==0 && ~isempty(e)
   Ds = ltipack.ssdata(a,b,c,d,e,0);  % D may contain residual internal delays
   Ds.Scaled = true;
   [isP,Ds] = isproper(Ds);
   if isP
      % System is proper. Use reduced model with nonsingular E
      [a,b,c,d,~,e] = getABCDE(Ds);
   elseif isempty(fBand) || isinf(fBand(2))
      % Gain is infinite at w=Inf for improper systems
      gpeak = Inf;  fpeak = Inf;  return
   end
end

% Compute norm
if allfinite(d)
   % Compute equivalent continuous-time poles and test frequencies
   if isempty(e)
      r = eig(a);
   else
      r = eig(a,e);  % E may be singular in DT
   end
   if Ts>0
      r = log(r)/Ts;
   end
   fTest = abs(r(isfinite(r)));
   fTest = unique([-fTest ; 0 ; fTest]);  % for complex data
   [gpeak,fpeak] = ltipack.norminf(a,b,c,d,e,Ts,tol,Inf,fTest,fBand);
else
   gpeak = NaN;  fpeak = NaN;
end

