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

%   Author(s): P. Gahinet
%   Copyright 1986-2000 The MathWorks, Inc. 
%   $Revision: 1.5 $ $Date: 2000/06/14 20:09:11 $

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

% 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;

% Handle different modes
switch lower(Mode)
case {'tight' 'alltight'}
   % 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.
   VisFilter = {};
   if strcmp(lower(Mode),'tight'),
      VisFilter = {'Visible','on'};
   end
   
   % Determine data extent for each axis
   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)] = UpdateExtent1(Xlims(:,k),Ylims(:,k),...
            get(ResponseLines(i),'Xdata'),get(ResponseLines(i),'Ydata'));
      end
   end
   
   % Derive overall axis limits
   Xlims = GetLims('X',min(Xlims(1,:,:),[],2),max(Xlims(2,:,:),[],2),AllRespProps.Xlim);
   Ymax = max(Ylims(2,:,:),[],3);
   Ylims = GetLims('Y',-Ymax,Ymax,AllRespProps.Ylim);
     
case 'critical'
   % Zoom in around the critical point
   MinExtent = 4;  % +/-4 in each direction away from -1,0 for Mode=CriticalPoint
   MinRadius = zeros(ny,nu);
   jay = sqrt(-1);
   
   % For each plot axis, determine the smallest ball centered at (-1,0)
   % containing data for all visible models
   for k=1:ny*nu,
      ResponseLines = findobj(get(PlotAxes(k),'Children'),...
                              'Tag','LTIresponseLines','Visible','on');
      for i=1:length(ResponseLines),
         freqdata = get(ResponseLines(i),'Xdata') + jay * get(ResponseLines(i),'Ydata');
         if any(diff(freqdata)),  % ignore static gains
            mingap = min(abs(1+freqdata));
            MinRadius(k) = max(MinRadius(k),mingap);
         end
      end
   end
   
   % Infer desirable X and Y ranges
   Xextents = max(MinExtent,1.5*max(MinRadius,[],1));  % 1xnu
   Yextents = max(MinExtent,1.5*max(MinRadius,[],2));  % nyx1
   Xbounds = num2cell([-1-Xextents(:) , -1+Xextents(:)],2);
   Ybounds = num2cell([-Yextents(:) , Yextents(:)],2);
   
   % Adjust limits to make sure they don't exceed overall data extent
   for j=1:nu,
      for i=1:ny,
         ResponseLines = findobj(get(PlotAxes(i,j),'children'),...
            'Tag','LTIresponseLines','Visible','on');
         % Determine line data extent within ranges XBOUNDS and YBOUNDS
         for k=1:length(ResponseLines)
            [Xlims(:,i,j),Ylims(:,i,j)] = UpdateExtent2(Xlims(:,i,j),Ylims(:,i,j),...
               get(ResponseLines(k),'Xdata'),get(ResponseLines(k),'Ydata'),...
               Xbounds{j},Ybounds{i});
         end
      end
   end
   
   % Set X limits
   Xlims = GetLims('X',min(Xlims(1,:,:),[],2),max(Xlims(2,:,:),[],2),Xbounds);
   for j=1:nu,
      Xlims{j} = [max(Xbounds{j}(1),Xlims{j}(1)) , min(Xbounds{j}(2),Xlims{j}(2))];
   end
   
   % Set Y limits
   Ylims = GetLims('Y',min(Ylims(1,:,:),[],3),max(Ylims(2,:,:),[],3),Ybounds);
   for i=1:ny,
      Ylims{i}(1) = max(Ybounds{i}(1),Ylims{i}(1));
      Ylims{i}(2) = -Ylims{i}(1);
   end      
end


%==========================
% Subfunction UpdateExtent1
%==========================
function [Xlims,Ylims] = UpdateExtent1(Xlims,Ylims,Xdata,Ydata)
%   Updates extent XLIMS and YLIMS.  Applies clipping algorithm to eliminate
%   large branches associated with poles near the imaginary axis or unit circle.

fin = isfinite(Xdata) & isfinite(Ydata);
Xdata = Xdata(fin);
Ydata = Ydata(fin);
Xmin = min(Xdata);
Xmax = max(Xdata);
Ymin = min(Ydata);
Ymax = max(Ydata);
LCP = any(Xdata<-1 & Ydata==0);
RCP = any(Xdata>-1 & Ydata==0);

% Clip unbounded branches associated with poles near the imaginary axis
% or unit circle
if Xmax-Xmin>20 | Ymax-Ymin>20,
   % Compute gain and phase of frequency response 1+H
   H = Xdata+j*Ydata;
   G = abs(1+H);
   P = unwrap(angle(1+H));
   % Find peak values of |1+H|
   dG = diff(G);
   Peak1 = (dG(1:end-1)>0 & dG(2:end)<0);  % single-point peak
   Peak2 = (dG(1:end-2)>0 & dG(3:end)<0);  % double-point peak (jw axis mode)
   ipeak = find(G>50 & [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 = 1+min(G);
   for i=1:npeak,
      Discard{i} = find(abs(P-P(ipeak(i)))<=pi/90 & G>10*Gmin);
   end
   Discard = cat(2,Discard{:});
   Xdata(:,Discard) = [];
   Ydata(:,Discard) = [];
   % Compute reduced bounding box
   Xminr = min(Xdata);  Xmaxr = max(Xdata);  Xextent = max(1,Xmaxr-Xminr);
   Yminr = min(Ydata);  Ymaxr = max(Ydata);  Yextent = max(1,Ymaxr-Yminr);
   Xmin = max(Xmin,Xminr-0.25*Xextent);
   Xmax = min(Xmax,Xmaxr+0.25*Xextent);
   Ymin = max(Ymin,Yminr-0.25*Yextent);
   Ymax = min(Ymax,Ymaxr+0.25*Yextent);
end

% Make sure Xmin<-1 when there's real data left of (-1,0) (same for Xmax)
if LCP,
   Xmin = min(Xmin,-2);
end
if RCP,
   Xmax = max(Xmax,0);
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 UpdateExtent2
%==========================
function [Xlims,Ylims] = UpdateExtent2(Xlims,Ylims,Xdata,Ydata,Xbounds,Ybounds)
%   Updates extent XLIMS and YLIMS by incorporating data within range
%   XBOUNDS and YBOUNDS.

% Exclude data outside of X-Y box specified by XBOUNDS and YBOUNDS
% Keep extra point to left and right of excluded portions to guard against 
% sparsely distributed points along unbounded branches. 
% Cf. sys = (tf(100,[1 0.001 100]) + tf(10,[1 1 10])) * tf(1,[1 1],'iod',.1);
InScope = find(Xdata>=Xbounds(1) & Xdata<=Xbounds(2) & Ydata>=Ybounds(1) & Ydata<=Ybounds(2));
InScope = sort([InScope , InScope(diff([1 InScope])>1)-1 ,...
                          InScope(diff([InScope length(Xdata)])>1)+1]);

% Update current X and Y limits
Xlims(1) = min([Xlims(1),min(Xdata(InScope))]);
Xlims(2) = max([Xlims(2),max(Xdata(InScope))]);
Ylims(1) = min([Ylims(1),min(Ydata(InScope))]);
Ylims(2) = max([Ylims(2),max(Ydata(InScope))]);

% Make sure Xmin<-1 when there's real data left of (-1,0) (same for Xmax)
if any(Xdata<-1 & Ydata==0),
   Xlims(1) = min(Xlims(1),-2);
end
if any(Xdata>-1 & Ydata==0),
   Xlims(2) = max(Xlims(2),0);
end


%====================
% Subfunction GetLims
%====================
function Lims = GetLims(XY,Min,Max,CurrentLims)
%GETLIMS  Polishes computed axis limits.

CP = -strcmpi(XY,'x');  % crtical point coordinate

% Set undetermined X limits to current limits
for j=find(Min>Max),
   Min(j) = CurrentLims{j}(1);
   Max(j) = CurrentLims{j}(2);
end

% Ensure minimum extent of 0.2, and slightly increase current extent
Median = (Min+Max)/2;
Extent = Max-min(CP,Min);
lcp = (Min<=CP);
Min(~lcp) = min(CP,Median(~lcp)-0.1);  % include critical point
Min(lcp) = min(Min(lcp)-0.03*Extent(lcp),Median(lcp)-0.1);
Max = max(Max+0.03*Extent,Median+0.1);
   
% Format outputs, include critical point
Lims = num2cell([Min(:) Max(:)],2);



      