function [gridlines,labels] = zpchart(varargin)
%ZPCHART  Generates the z-plane grid lines
%
%    [GRIDLINES,LABELS] = ZPCHART(AX) plots z-plane grid lines 
%    on the axes AX.  The range and spacing of the natural 
%    frequencies and damping ratios is determined automatically
%    based on the axes limits of AX.  GRIDLINES and LABELS
%    contain the handles for the lines and labels of the plotted
%    grid.
%
%    [GRIDLINES,LABELS] = ZPCHART(AX,ZETA,WN) plots an z-plane 
%    grid on the axes AX for the set of damping ratios and natural 
%    frequencies specified in ZETA and WN, respectively.
%
%    Additional grid options may be provided in a structure as the
%    final input argument, as follows:
%
%       [GRIDLINES,LABELS] = ZPCHART(AX,s)
%    or
%       [GRIDLINES,LABELS] = ZPCHART(AX,ZETA,WN,s)
%
%    where s is a structure with the following fields:
%
%       s.FrequencyUnits = '[ Hz | {rad/sec} ]';
%       s.GridLabelType  = '[ {damping} | overshoot ]';
%       s.SampleTime     = '[ real scalar | {-1} ]';
%
%    Note that if frequency units of 'Hz' have been specified and
%    WN have been provided, then it is assumed that the WN values
%    are provided in units of 'Hz'.
%
%    See also PZMAP, ZGRID.

%   Revised: Adam DiVergilio, 11-99
%   Copyright 1986-2000 The MathWorks, Inc. 
%   $Revision: 1.16 $  $Date: 2000/09/07 18:54:49 $

%---Initialize return values
gridlines = [];
labels = [];

%---If final argument is a structure, then use options provided
if isstruct(varargin{end})
   s = varargin{end};
   FrequencyUnits = s.FrequencyUnits;
   GridLabelType  = s.GridLabelType;
   SampleTime     = s.SampleTime;
else
   FrequencyUnits = 'rad/sec';
   GridLabelType  = 'damping';
   SampleTime     = -1;
end

%---Axes info
ax = varargin{1};
XLIM = get(ax,'XLim');
YLIM = get(ax,'YLim');
dXLIM = diff(XLIM);
dYLIM = diff(YLIM);
ContainsXAxis = YLIM(1)<=0 & YLIM(2)>=0;
ContainsYAxis = XLIM(1)<=0 & XLIM(2)>=0;
ContainsOrigin = ContainsXAxis & ContainsYAxis;
Color = get(ax,'XColor');
LColor = Color;
TColor = localTextColor(ax);
FontSize = get(ax,'FontSize');
FontWeight = get(ax,'FontWeight');
FontAngle = get(ax,'FontAngle');

%---Plot unit circle (always)
t = 0:2*pi/100:2*pi;
unitcircle = line('Parent',ax,...
   'XData',sin(t),'YData',cos(t),...
   'LineStyle','-','Color',LColor,...
   'Tag','CSTgridLines','HitTest','off','HandleVisibility','off','HelpTopicKey','isodampinggrid');

%---User-defined damping ratios and natural frequencies
if nargin>1 & isnumeric(varargin{2})
   zeta = varargin{2};
   wn = unitconv(varargin{3},FrequencyUnits,'rad/sec');  % Work in 'rad/sec'
   %---Make sure we are working with sorted rows
   zeta = sort(zeta(:)');
   wn = sort(wn(:)');
   %---Drop bogus values
   zeta = zeta(:,(zeta>=0 & zeta<=1));
   wn = wn(:,(wn>0 & wn<=pi));
   %---Round extremely small zeta values to zero
   zeta(:,zeta<10*eps) = 0;
   %---Must have at least one zeta and one wn value to continue
   if isempty(zeta) | isempty(wn)
      %---Return handle of unit circle only
      gridlines = unitcircle(:);
      return
   end
   %---Build strings for wn values (used in place of pi/T fractions)
   wns = cell(length(wn),1);
   for n=1:length(wn)
      wns{n} = sprintf('%0.3g',unitconv(wn(n),'rad/sec',FrequencyUnits));  % Show in requested units
   end

%---Default damping ratios and natural frequencies
else
   zeta = .1:.1:.9;
   wn = pi/10:pi/10:pi;
   if SampleTime>0
      for n=1:length(wn)
         wns{n} = sprintf('%0.3g',unitconv(wn(n)/SampleTime,'rad/sec',FrequencyUnits));
      end
   else
      if strcmpi(FrequencyUnits,'Hz')
         wns = {'0.05/T','0.10/T','0.15/T','0.20/T','0.25/T',
                '0.30/T','0.35/T','0.40/T','0.45/T','0.50/T'};
      else
         wns = {'\pi/10T','\pi/5T','3\pi/10T','2\pi/5T','\pi/2T',
                '3\pi/5T','7\pi/10T','4\pi/5T','9\pi/10T','\pi/T'};
      end
   end
end
wns = strrep(wns,'+00','');
wns = strrep(wns,'+0','');

%---Plot zeta lines (damping)
m = tan(asin(zeta)) + i;
zzeta = [exp((0:pi/100:pi)'*(-m)); NaN*ones(size(zeta))];
zzeta = zzeta(:);
rzzeta = real(zzeta);
izzeta = imag(zzeta);
zetalines = line('Parent',ax,...
   'XData',[rzzeta; rzzeta],'YData',[izzeta; -izzeta],...
   'LineStyle',':','Color',LColor,...
   'Tag','CSTgridLines','HitTest','off','HandleVisibility','off','HelpTopicKey','isodampinggrid');

%---Plot wn lines (natural frequency)
e_itheta = exp(i * (pi/2:pi/50:pi)');
e_r = exp(wn);
zwn = [(ones(length(e_itheta),1)*e_r).^(e_itheta*ones(size(e_r))); NaN*ones(size(e_r))];
zwnc = zwn(:);
rzwn = real(zwnc);
izwn = imag(zwnc);
wnlines = line('Parent',ax,...
   'XData',[rzwn; rzwn],'YData',[izwn; -izwn],...
   'LineStyle',':','Color',LColor,...
   'Tag','CSTgridLines','HitTest','off','HandleVisibility','off','HelpTopicKey','isofrequencygrid');

%---Return handles of new gridlines
gridlines = [wnlines(:); zetalines(:); unitcircle(:)];

%---Limits may have changed due to axis autoscaling,
%--- so requery limits prior to placing text
XLIM = get(ax,'XLim');
YLIM = get(ax,'YLim');
dXLIM = diff(XLIM);
dYLIM = diff(YLIM);
ContainsXAxis = YLIM(1)<=0 & YLIM(2)>=0;
ContainsYAxis = XLIM(1)<=0 & XLIM(2)>=0;
ContainsOrigin = ContainsXAxis & ContainsYAxis;

%---Construct wn midpoint information (for zeta label positioning)
%
%---Create wn vector which contains 0,pi
if wn(end)==pi
   wntmp = [0 wn];
else
   wntmp = [0 wn pi];
end
%---Vector of midpoints between plotted wn lines
wnmid = (wntmp(1:end-1)+wntmp(2:end))/2;
%---Coordinates of wnmid on z-plane
e_itheta = exp(i * (pi/2:pi/50:pi)');
e_r = exp(wnmid);
zwnmid = [(ones(length(e_itheta),1)*e_r).^(e_itheta*ones(size(e_r)))];

%---If y=0 line is showing, draw frequency labels above and below x-axis
if ContainsXAxis
   wn = [wn wn];
   zwn = [zwn conj(zwn)];
   wns = [wns wns];
   wnmid = [wnmid wnmid];
   zwnmid = [zwnmid conj(zwnmid)];
%---If only negative y-values are showing, only draw frequency labels below x-axis
elseif YLIM(2)<0
   zwn = conj(zwn);
   zwnmid = conj(zwnmid);
end

%---Space fraction from axes edge
edge = 0.015;

%---Space fraction of axes used to determine label crowding
crowdlim = 1/6;

%---Try to find a label position for each wn value
for n=1:length(wn)
   %---Find coordinates of wn line which lie within axes limits
   wnz = zwn(~isnan(zwn(:,n)),n);
   idx = find((real(wnz)>=XLIM(1))&(real(wnz)<=XLIM(2))&...
              (imag(wnz)>=YLIM(1))&(imag(wnz)<=YLIM(2)));
   %---Initialize label position
   wnx = [];
   wny = [];
   %---If this wn line is visible, add label at position set by wn and axes limits
   if ~isempty(idx)
      %---Visible coordinate along wn which is closest to unit circle
      firstvis = idx(1);
      %---If unit circle is shown at this wn, use it for label position
      if firstvis==1
         wnx = real(wnz(1));
         wny = imag(wnz(1));
         dx = real(wnz(2))-wnx;
         dy = imag(wnz(2))-wny;
         slp = dy/dx;
         if     wn(n)>pi/2, HA = 'left';
         elseif wn(n)<pi/2, HA = 'right';
         else,              HA = 'center';
         end
      %---Otherwise, interpolate a point along axes edge
      else
         w1 = wnz(firstvis-1);
         w2 = wnz(firstvis);
         x = [real(w1) real(w2)];
         y = [imag(w1) imag(w2)];
         dx = diff(x);
         dy = diff(y);
         slp = dy/dx;
         if y(1)>(YLIM(2)-dYLIM*edge)
            if ~(ContainsXAxis & abs(YLIM(2)/dYLIM)<crowdlim)
               wny = YLIM(2)-dYLIM*edge;
               wnx = x(1) + (wny-y(1))/slp;
               HA = 'center';
            end
         elseif y(1)<(YLIM(1)+dYLIM*edge)
            if ~(ContainsXAxis & abs(YLIM(1)/dYLIM)<crowdlim)
               wny = YLIM(1)+dYLIM*edge;
               wnx = x(1) + (wny-y(1))/slp;
               HA = 'center';
            end
         elseif x(1)<XLIM(1)+dXLIM*edge;
            if ~(ContainsYAxis & abs(XLIM(1)/dXLIM)<crowdlim)
               wnx = XLIM(1)+dXLIM*edge;
               wny = y(1) + (wnx-x(1))*slp;
               HA = 'left';
            end
         else %--- x(1)>XLIM(2)-dXLIM*edge
            if ~(ContainsYAxis & abs(XLIM(2)/dXLIM)<crowdlim)
               wnx = XLIM(2)-dXLIM*edge;
               wny = y(1) + (wnx-x(1))*slp;
               HA = 'right';
            end
         end
      end
      %---Select vertical alignment based on dy (slope in y)
      if dy<=0, VA = 'top';
      else,     VA = 'bottom';
      end
   end
   %---Draw label if valid coordinate was determined
   if ~isempty(wnx)
      labels(length(labels)+1) = text(...
         'Units','data','Position',[wnx wny],'String',wns(n),...
         'Parent',ax,'Color',TColor,'FontSize',FontSize,'FontWeight',FontWeight,'FontAngle',FontAngle,...
         'Clipping','on','Rotation',0,'HorizontalAlignment',HA,'VerticalAlignment',VA,...
         'Tag','CSTgridLines','HitTest','off','HandleVisibility','off','HelpTopicKey','isofrequencygrid');
   end
end

%---If the default locations for zeta labels is visible, use them
wndef = 3.5*pi/10;
zzdef = exp(wndef*(-zeta + i*sqrt(1-zeta.^2)));
zzdefr = real(zzdef);
zzdefi = imag(zzdef);
zzdefin = -zzdefi;
if all(zzdefr>=XLIM(1) & zzdefr<=XLIM(2) & zzdefi>=YLIM(1) & zzdefi<=YLIM(2)) | ...
   all(zzdefr>=XLIM(1) & zzdefr<=XLIM(2) & zzdefin>=YLIM(1) & zzdefin<=YLIM(2))
   wnzeta = wndef;

%---Otherwise, select most-visible vector from wnmid values
else
   %---Create vector of values representing visible zeta range of each wnmid vector
   for n=1:length(wnmid)
      %---Find coordinates of wnmid line which lie within axes limits
      wnz = zwnmid(:,n);
      idx = find((real(wnz)>=XLIM(1))&(real(wnz)<=XLIM(2))&...
                 (imag(wnz)>=YLIM(1))&(imag(wnz)<=YLIM(2)));
      if isempty(idx)
         zdiff(n) = 0;
      else
         %---Visible coordinate along wnmid which is closest to unit circle
         zmin = real((-((log(wnz(idx(1)))).^2)-(wnmid(n)).^2)./(2*wnmid(n)*log(wnz(idx(1)))));
         %---Visible coordinate along wnmid which is farthest from unit circle
         zmax = real((-((log(wnz(idx(end)))).^2)-(wnmid(n)).^2)./(2*wnmid(n)*log(wnz(idx(end)))));
         %---Fix round-off error
         if abs(zmin)<10*eps,   zmin = 0; end
         if abs(1-zmax)<10*eps, zmax = 1; end
         %---Visible zeta range
         zdiff(n) = zmax-zmin;
      end
   end
   %---wnmid value associated with max(zdiff) is where the zeta labels should go
   [tmp,zidx] = max(zdiff);
   wnzeta = wnmid(zidx);
end

%---Draw zeta labels
if abs(YLIM(2))>=abs(YLIM(1))
   wnzetasign = 1;
else
   wnzetasign = -1;
end
zz = exp(wnzeta*(-zeta + i*sqrt(1-zeta.^2)));
if strcmpi(GridLabelType,'overshoot')
   dispzeta = exp(-(zeta*pi)./sqrt(1-zeta.*zeta));
   dispzeta = round(1e5*dispzeta)/1e3; % round off small values
else
   dispzeta = zeta;
end
for n=1:length(zeta)
   labels(length(labels)+1) = text(...
      'Units','data','Position',[real(zz(n)) wnzetasign*imag(zz(n))],...
      'String',sprintf('%0.3g',dispzeta(n)),...
      'Parent',ax,'Color',TColor,'FontSize',FontSize,'FontWeight',FontWeight,'FontAngle',FontAngle,...
      'Clipping','on','Rotation',0,'HorizontalAlignment','center','VerticalAlignment','middle',...
      'Tag','CSTgridLines','HitTest','off','HandleVisibility','off','HelpTopicKey','isodampinggrid');
end

%---Return handles of new labels
labels = labels(:);

%---Restack axes children so that grid lines/labels are on bottom
AllKids = allchild(ax);
BotKids = [labels;gridlines];
TopKids = AllKids(~strcmp(get(AllKids,'Tag'),'CSTgridLines'));
set(ax,'Children',[TopKids(:);BotKids(:)]);


%%%%%%%%%%%%%%%%%%
% localTextColor %
%%%%%%%%%%%%%%%%%%
function C = localTextColor(Ax)
%---Determine good color for grid text in axes Ax
 C  = get(Ax,'XColor');
 AC = get(Ax,'Color');
 if isstr(AC), AC = get(get(Ax,'Parent'),'Color'); end
 if sum(AC)>1.5
    C = C*0.5;
 else
    C = C + 0.5*([1 1 1]-C);
 end
