function h0 = nicholsplot(varargin)
%NICHOLSPLOT  Nichols frequency response of linear systems.
%
%   NICHOLSPLOT, an extension of NICHOLS, provides a command line interface
%   for customizing the plot appearance.
%
%   NICHOLSPLOT(SYS) draws the Nichols plot of the dynamic system SYS.  The
%   frequency range and number of points are chosen automatically.  See
%   BODE for details on the notion of frequency in discrete-time.
%
%   NICHOLSPLOT(SYS,{WMIN,WMAX}) draws the Nichols plot for frequencies
%   between WMIN and WMAX (in radians/timeunit, where timeunit is specified
%   in SYS.TimeUnit).
%
%   NICHOLSPLOT(SYS,W) uses the user-supplied vector W of frequencies (in
%   radians/timeunit, where timeunit is specified in SYS.TimeUnit) at which
%   the Nichols response is to be evaluated. See LOGSPACE to generate
%   logarithmically spaced frequency vectors.
%
%   NICHOLSPLOT(SYS1,SYS2,...,W) plots the Nichols plot of multiple LTI
%   models SYS1,SYS2,... on a single plot.  The frequency vector W
%   is optional.  You can also specify a color, line style, and marker 
%   for each system, as in  
%      nicholsplot(sys1,'r',sys2,'y--',sys3,'gx').
%
%   NICHOLSPLOT(AX,...) plots into the axes with handle AX.
%
%   NICHOLSPLOT(..., PLOTOPTIONS) plots the Nichols chart with the options
%   specified in PLOTOPTIONS. See NICHOLSOPTIONS for more details. 
%
%   H = NICHOLSPLOT(...) returns the handle to the Nichols plot. You can 
%   use this handle to customize the plot with the GETOPTIONS and 
%   SETOPTIONS commands.  See NICHOLSOPTIONS for a list of available plot 
%   options.
%
%   Example:
%       sys = rss(5);
%       h = nicholsplot(sys);
%       % Change units to Hz 
%       setoptions(h,'FreqUnits','Hz');
%
%   See also NICHOLS, NICHOLSOPTIONS, WRFC/SETOPTIONS, WRFC/GETOPTIONS, DYNAMICSYSTEM.

%   Authors: P. Gahinet, B. Eryilmaz
%   Copyright 1986-2010 The MathWorks, Inc.

% Get argument names
for ct = length(varargin):-1:1
    ArgNames(ct,1) = {inputname(ct)};
end

% Parse input list
% Check for axes argument
if ishghandle(varargin{1},'axes')
    ax = varargin{1};
    varargin(1) = [];
    ArgNames(1) = [];
else
    ax = [];
end

try
    [sysList,Extras,OptionsObject] = DynamicSystem.parseRespFcnInputs(varargin,ArgNames);
    if ~isempty(OptionsObject) && ~isa(OptionsObject,'plotopts.NicholsPlotOptions')
       error('Controllib:plots:InvalidPlotOptions',...
                getString(message('Controllib:plots:InvalidPlotOptions','nicholsoptions')));
   end
    [sysList,w] = DynamicSystem.checkBodeInputs(sysList,Extras);
    TimeUnits = sysList(1).System.TimeUnit; % first system determines units
    % Derive plot I/O size and I/O names
    [InputName,OutputName] = mrgios(sysList.System);
    % Warn about and skip empty systems
    isEmptySys = arrayfun(@(x) isempty(x.System),sysList);
    if any(isEmptySys)
        ctrlMsgUtils.warning('Control:analysis:PlotEmptyModel')
    end
    sysList = sysList(~isEmptySys);
    % Check time unit consistency when specifying w or {wmin,wmax}
    if ~(isempty(w) || ltipack.hasMatchingTimeUnits(TimeUnits,sysList.System))
        ctrlMsgUtils.error('Control:analysis:AmbiguousFreqSpec')
    end
    
    % Create LTISource and check compatibility
    nsys = numel(sysList);
    src = repmat(handle(NaN),nsys,1);
    style = cell(nsys,1);
    for ct=1:length(sysList)
        sysInfo = sysList(ct);
        src(ct) = resppack.ltisource(sysInfo.System, 'Name', sysInfo.Name);
        checkCompatibility(src(ct),'nichols');
        style{ct} = sysInfo.Style;
    end
    
    % Nichols plot
    % Create plot (visibility ='off')
    if isempty(ax)
        ax = gca;
    end
    h = ltiplot(ax,'nichols',InputName,OutputName,OptionsObject,cstprefs.tbxprefs);
catch ME
    throw(ME)
end

if isnumeric(w)
    % Unique frequencies, to avoid interpolation incompatibility
    % for other calculations. Resolution for G154921
    w = unique(w);
end

% Configure responses
for ct=1:length(src)
      % Link each response to system source
      r = h.addresponse(src(ct));
      DefinedCharacteristics = src(ct).getCharacteristics('nichols');
      r.setCharacteristics(DefinedCharacteristics);
      r.DataFcn = {'magphaseresp' src(ct) 'nichols' r w};
      % Styles and preferences
      initsysresp(r,'nichols',h.Options,style{ct})
end


% Draw now
if strcmp(h.AxesGrid.NextPlot, 'replace')
    h.Visible = 'on';  % new plot created with Visible='off'
else
    draw(h)  % hold mode
end

% Right-click menus
ltiplotmenu(h, 'nichols');

if nargout
    h0 = h;
end
