function [Xlims,Ylims] = niclims(NicRespObj,Mode)
%NICLIMS  Calculates axis limits for Nyquist plots.
%
%   [XLIMS,YLIMS] = NICLIMS(NICRESPOBJ,MODE) computes desirable axis
%   limits for Nyquist plots.  NICRESPOBJ is a Nichols Response Object,
%   and MODE selects between
%      'tight'   :  Tight axis limits for visible models
%      'alltight':  Same but including all models (visible or not)
%   The outputs XLIMS and YLIMS are NUx1 and NYx1 cell arrays of
%   X and Y limits.
%
%   See also AXESLIMS, PLOTGRID.

%   Authors: P. Gahinet
%   Revised: A. DiVergilio
%   Copyright 1986-2000 The MathWorks, Inc. 
%   $Revision: 1.7 $ $Date: 2000/06/14 20:09:39 $

% Get plot axes
AllRespProps = get(NicRespObj,'Undocumented');
PlotAxes = AllRespProps.PlotAxes;
[ny,nu] = size(PlotAxes);

% Unit/Scale info
ISLOGX = 0;
ISLOGY = strcmpi(getscale(NicRespObj,'MagnitudeScale'),'log');
DB     = strcmpi(AllRespProps.MagnitudeUnits,'dB');
DEG    = strcmpi(AllRespProps.PhaseUnits,'deg');

% Initialize X and Y ranges
Xlims = zeros(2,ny,nu);  Xlims(1,:) = Inf;  Xlims(2,:) = -Inf;
Ylims = zeros(2,ny,nu);  Ylims(1,:) = Inf;  Ylims(2,:) = -Inf;

% Set include switch (all vs. visible only)
VisFilter = {};
if strcmp(lower(Mode),'tight')
   VisFilter = {'Visible','on'};
end

% Determine axis limits for visible models (zoom option 'tight')
% or entire set of models (used for grid computation).
% RE: XlimMode and YlimMode are ignored in both cases.
for k=1:ny*nu
   % Find response lines
   ResponseLines = findobj(get(PlotAxes(k),'Children'),'Tag','LTIresponseLines',VisFilter{:});
   % Determine line data extent. Clip infinite branches due to poles near imag. axis
   for i=1:length(ResponseLines)
      [Xlims(:,k),Ylims(:,k)] = UpdateExtent(Xlims(:,k),Ylims(:,k),...
         get(ResponseLines(i),'XData'),get(ResponseLines(i),'YData'),DB,DEG,ISLOGY);
   end
end

%---Derive overall X-axis limits
Xmin = min(Xlims(1,:,:),[],2);
Xmax = max(Xlims(2,:,:),[],2);
Xlims = num2cell([Xmin(:) Xmax(:)],2);
if DEG
   Step = 30;
   MinExtent = 2;
   MinValue = -Inf;
else
   Step = 0.5;
   MinExtent = .05;
   MinValue = -Inf;
end
for k=1:nu
   Xlims{k} = PolishLims(Xlims{k},AllRespProps.Xlim{k},Step,MinExtent,ISLOGX,MinValue);
end

%---Derive overall Y-axis limits
Ymin = min(Ylims(1,:,:),[],3);
Ymax = max(Ylims(2,:,:),[],3);
Ylims = num2cell([Ymin(:) Ymax(:)],2);
if DB
   Step = 10;
   MinExtent = 2;
   MinValue = -Inf;
else
   Step = 0.5;
   MinExtent = 0.5;
   if ISLOGY
      MinValue = eps;
   else
      MinValue = 0;
   end
end
for k=1:ny
   Ylims{k} = PolishLims(Ylims{k},AllRespProps.Ylim{k},Step,MinExtent,ISLOGY,MinValue);
end


%=========================
% Subfunction UpdateExtent
%=========================
function [Xlims,Ylims] = UpdateExtent(Xlims,Ylims,Xdata,Ydata,DB,DEG,ISLOGY)
%   Updates extent XLIMS and YLIMS by incorporating data within range
%   XBOUNDS and YBOUNDS.  Applies clipping algorithm to eliminate
%   large branches associated with poles near the imaginary axis
%   or unit circle.

if DEG
   Xdel = 10;
   Xsmall = 2;
else
   Xdel = 0.2;
   Xsmall = 0.05;
end
if DB
   Ydel = 10;
   Ybig = 40;
   Ybig1 = -40;
   Ybig2 = 40;
   Y25 = 25;
else
   Ydel = 3;
   Ybig = 100;
   Ybig1 = .01;
   Ybig2 = 100;
   Y25 = unitconv(25,'dB','abs');
end

fin = isfinite(Xdata) & isfinite(Ydata);
Xdata = Xdata(fin);
Ydata = Ydata(fin);
Xmin = min(Xdata);
Xmax = max(Xdata);
Ymin = min(Ydata);
Ymax = max(Ydata);

% Clip unbounded branches associated with poles near the imaginary axis
% or unit circle
if Ymax-Ymin>Ybig & (Ymin<Ybig1 | Ymax>Ybig2)
   % Find peak gain values
   Gabs = abs(Ydata);
   dG = diff(Ydata);
   Peak1 = (dG(1:end-1)>0 & dG(2:end)<0);   % single-point peak
   Peak2 = (dG(1:end-2)>0 & dG(2:end-1)<0); % double-point peak (jw axis mode)
   ipeak = find(Gabs>Ybig & [1 , (Peak1 | [Peak2 0] | [0 Peak2]) , 1]);
   % Clip unbounded branches by deleting data with phase within 2 degrees
   % of phase at peak gain values
   npeak = length(ipeak);
   Discard = cell(npeak,1);
   Gmin = min(Gabs);
   for i=1:npeak
      Discard{i} = find(abs(Xdata-Xdata(ipeak(i)))<=Xsmall & Gabs>=Y25+Gmin);
   end
   Discard = cat(2,Discard{:});
   Xdata(:,Discard) = [];
   Ydata(:,Discard) = [];
   % Compute reduced bounding box
   Xminr = min(Xdata);   Xmaxr = max(Xdata);
   Yminr = min(Ydata);   Ymaxr = max(Ydata);
   Xmin = max(Xmin,Xminr-Xdel);
   Xmax = min(Xmax,Xmaxr+Xdel);
   Ymin = max(Ymin,Yminr-Ydel);
   Ymax = min(Ymax,Ymaxr+Ydel);
end

% Update current X and Y limits
Xlims(1) = min([Xlims(1),Xmin]);
Xlims(2) = max([Xlims(2),Xmax]);
Ylims(1) = min([Ylims(1),Ymin]);
Ylims(2) = max([Ylims(2),Ymax]);


%=======================
% Subfunction PolishLims
%=======================
function Lims = PolishLims(Lims,CurrentLims,Step,MinExtent,ISLOG,MinValue)
%POLISHLIMS  Polishes computed axis limits.

%---Set undetermined limits to current limits
if Lims(1)>Lims(2)
   Lims = CurrentLims;
end

%---Save lower limit incase we need to adjust for log scale
LowerLim = Lims(1);

%---Ensure minimum extent is met
if (Lims(2)-Lims(1))<MinExtent
   Median = (Lims(1)+Lims(2))/2;
   Lims = Median + [-1 1]*MinExtent/2;
end

%---Round out to nearest Step locations
if (Lims(2)-Lims(1))>Step
   Lims(1) = Step * floor(Lims(1)/Step);
   Lims(2) = Step * ceil(Lims(2)/Step);
else
   Lims(1) = MinExtent * floor(Lims(1)/MinExtent);
   Lims(2) = MinExtent * ceil(Lims(2)/MinExtent);
end

%---If scale is log and lower limit is not >=0, then adjust lower limit
if ISLOG & Lims(1)<=0
   Lims(1) = LowerLim;
end

%---Make sure limit lower limit is >=MinValue
Lims(1) = max(MinValue,Lims(1));
