function [varargout] = hsvplot(varargin)
%HSVPLOT  Plots the Hankel singular values of an LTI model.
%
%   HSVPLOT(SYS) plots the Hankel singular values (HSVs) of the LTI model 
%   SYS. See BALRED for details on the meaning and purpose of HSVs. The
%   HSVs for the stable and unstable modes of SYS are shown in blue and 
%   red, respectively. Reduced-order models of various orders can be
%   obtained by dropping the states associated with the smallest HSVs.
%   HSVPLOT also shows the maximum approximation error as a function 
%   of the approximation order (number of states in reduced-order model).
%
%   HSVPLOT(AX,SYS,...) attaches the plot to the axes AX.
%
%   HSVPLOT(...,NUMOPTIONS,PLOTOPTIONS) specifies additional options for 
%   computing and plotting the results. Use BALREDOPTIONS to create 
%   NUMOPTIONS and HSVOPTIONS to create PLOTOPTIONS.
%
%   H = HSVPLOT(...) returns the handle H to the Hankel singular 
%   value plot. You can use this handle to customize the plot 
%   with the GETOPTIONS and SETOPTIONS commands. See HSVOPTIONS
%   for a list of available plot options.
%
%   Example:
%      sys = rss(20);
%      % Plot HSVs for relative-error method
%      numopt = balredOptions('ErrorBound','relative');
%      h = hsvplot(sys,numopt);
%      % Switch to linear scale
%      setoptions(h,'Yscale','linear')
%
%   See also HSVOPTIONS, BALREDOPTIONS, BALRED.

%	Author(s): P. Gahinet
%   Copyright 1986-2020 The MathWorks, Inc.
no = nargout;

% Check for axes argument
if nargin>0 && ishghandle(varargin{1})
   ax = varargin{1};
   varargin(1) = [];
else
   ax = -1;
end

% Get system
sys = varargin{1};

% Validate data and read options
if nmodels(sys)~=1
   error(message('Control:general:RequiresSingleModel','hsvplot'))
elseif any(iosize(sys)==0)
   % System without input or output
   error(message('Control:transformation:NotSupportedNoInputsorOutputs','hsvplot'))
end

% Options
ixn = find(cellfun(@(x) isa(x,'ltioptions.balred'),varargin));
ixp = find(cellfun(@(x) isa(x,'plotopts.HSVPlotOptions'),varargin));
if isempty(ixp)
   PlotOptions = hsvoptions();
   if isempty(ixn)
      % Watch for obsolete PV pair syntax:
      % hsvplot(sys,'offset',1.5,'AbsTol',1e-5,'RelTol',1e-4,'FreqIntervals',[0 2],'TimeIntervals',[3 5]);
      try
         NumOptions = balredOptions(varargin{2:end});
      catch ME
         error(ME.identifier,strrep(ME.message,'balred','hsvplot'))
      end
   else
      NumOptions = varargin{ixn};
   end
else
   PlotOptions = varargin{ixp};
   if isempty(ixn)
      % Backward compatibility: Remap BALRED options that used to live in hsvoptions
      NumOptions = balredOptions();
      NumOptions.AbsTol = PlotOptions.AbsTol;
      NumOptions.RelTol = PlotOptions.RelTol;
      NumOptions.Offset = PlotOptions.Offset;
      NumOptions.FreqIntervals = PlotOptions.FreqIntervals;
      NumOptions.TimeIntervals = PlotOptions.TimeIntervals;
   else
      NumOptions = varargin{ixn};
   end
end
NumOptions = validate(NumOptions);  % warn about incompatible options

% Compute HSV (needed to check computability)
try
   [~,Info] = balred(sys,NumOptions);
catch ME
   error(ME.identifier,strrep(ME.message,'balred','hsvplot'))
end

% Create plot (visibility ='off')
if ~ishghandle(ax,'axes')
   % Default to GCA
   % Call gca after argument parsing to prevent blank axis from appearing
   % if syntax error occurs
   ax = gca; %
end
try
   h = ltiplot(ax,'hsv',[],[],PlotOptions,cstprefs.tbxprefs);
catch ME
   throw(ME)
end

% Add HSV "response"
r = h.addresponse(1,1,1);
r.Data.HSV  = Info.HSV;
r.Data.ErrorBound  = [Info.ErrorBound;0];
r.Data.ErrorType = NumOptions.ErrorBound(1:3);

% Make plot visible
h.Visible = 'on';
legend(ax,'show')  % Legend on by default

% Add menu
h.addHSVMenu('yscale');
h.AxesGrid.addMenu('grid');
h.addMenu('fullview');

% Add properties menu
if usejava('MWT')
   h.addMenu('properties','Separator','on');
end

% Return handle if requested
if no>0
   varargout = {h};
end
