function BalData = hsvd(D,Options)
%HSVD  Computes low-rank approximations of gramians and corresponding
%      Hankel singular values.
%
%   Second input OPTIONS must be of class mor.SparseBalancedTruncationOptions.

%   Copyright 2023 The MathWorks, Inc.

% Note: CT only at this point
if hasInternalDelay(D)
   error(message('Control:transformation:BALROM14'))
elseif Options.Offset>0
   error(message('Control:transformation:BALROM30'))
end

% Rayleigh damping
if ~isempty(Options.Rayleigh)
   % Add implicit Rayleigh damping
   wn = Options.Rayleigh(1);
   zeta = Options.Rayleigh(2);
   alpha = zeta*wn;   beta = zeta/wn;
   nq = size(D.k,1);
   if isempty(D.m)
      D.c = D.c + alpha*speye(nq) + beta*D.k;
   else
      D.c = D.c + alpha*D.m + beta*D.k;
   end
elseif (ishermitian(D.m) && ishermitian(D.k) && all(diag(D.m)>=0) && all(diag(D.k)>=0))
   % Look for built-in Rayleigh damping and generate custom shifts when found
   RD = ltipack.util.RayleighDamping(D.c,D.m,D.k);
   if ~isempty(RD) && all(RD>0)
      % For shift generation
      Options.Rayleigh = [sqrt(RD(1)/RD(2)) sqrt(RD(1)*RD(2))];
   end
end

% Compute balancing data
BalData = hsvd(sparss(D),Options);
