function n = normh2(D)
% H2 norm of state-space model

%   Author(s): P. Gahinet
%   Copyright 1986-2018 The MathWorks, Inc.
Ts = D.Ts;

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

% NaN system
if ~allfinite(D.d)
   n = NaN;  return
elseif ~isstable(D)
   % Check stability before simplifying and scaling to avoid inconsistencies,
   % e.g., in H2SYN
   warning(message('Control:analysis:NormInfinite3'))
   n = Inf;  return
end

% Simplify and scale realization
% Note: 
%   1) In SYSTUNE context with loop openings, D may contain non-minimal 
%      poles at s=0 coming from, e.g., PIDs
%   2) Use getABCDE to properly handle I/O delays
D = sminreal(D);
if isreal(D.a) && isreal(D.b) && isreal(D.e)
   [a,b,c,d,~,e] = getABCDE(D);
else
   % Get explicit form (no complex descriptor version of Lyapunov solvers)
   [a,b,c,d] = getABCD(D);  e = [];
end
if ~D.Scaled
   [a,b,c,e] = ltipack.xscale(a,b,c,d,e,Ts,'Warn',false);
end

% Check properness and ensure E is invertible
if ~isempty(e)
   Ds = ltipack.ssdata(a,b,c,d,e,Ts);  % 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 Ts==0
      warning(message('Control:analysis:NormInfinite2'))
      n = Inf;
      return
   else
      error(message('Control:analysis:norm3'))
   end
end

% Compute norm
if Ts==0
   % Continuous H2 norm: 
   if norm(d,1)>0
      warning(message('Control:analysis:NormInfinite1'))
      n = Inf;
   else
      try
         % || G || = || R*c' ||_F where R'*R=P and a*P+P*a'+b*b' = 0
         R = lyapchol(a,b,e,'noscale');
         n = norm(R*c','fro');    % = sqrt(trace(c*P*c'))
      catch %#ok<*CTCH>
         warning(message('Control:analysis:NormInfinite3'))
         n = Inf;
      end
   end
else
   % Discrete H2 norm:
   try
      % || G ||^2 = trace(c*P*c'+d*d') where a*P*a'-P+b*b' = 0
      R = dlyapchol(a,b,e,'noscale');
      n = norm([R*c';d'],'fro');  % = sqrt(trace(c*P*c'+d*d'))
   catch
      warning(message('Control:analysis:NormInfinite3'))
      n = Inf;
   end
end
