function [phase,gain,labels] = nicchart(pmin,pmax,gmin)
%NICCHART  Generates the Nichols grid
%
%    [PHASE,GAIN,LABELS] = NICCHART(PMIN,PMAX,GMIN) generates the 
%    data for plotting the Nichols chart. PMIN and PMAX specify 
%    the phase interval in degrees and GMIN is the smallest gain 
%    in dB.  NICCHART returns
%       * PHASE and GAIN: grid data
%       * LABELS: three-row matrix containing the x,y positions 
%                 and the label values
%
%    See also NICHOLS, NGRID.

%	 Authors: K. Gondoly and P. Gahinet 
%	 Copyright 1986-2000 The MathWorks, Inc. 
%   $Revision: 1.3 $  $Date: 2000/06/14 19:59:11 $
%	 $Revision $  $Date: 2000/06/14 19:59:11 $

% Note: GMIN is assumed to be in dB  

j = sqrt(-1);
NanMat = NaN;

% Round Gmin from below to nearest multiple of -20dB,
% and Pmin,Pmax to nearest multiple of 360
gmin = min([-20 20*floor(gmin/20)]);
pmax = 360*ceil(pmax/360);
pmin = min(pmax-360,360*floor(pmin/360));

% (1) Generate isophase lines for following phase values:
p1 = [1 5 10 20 30 50 90 120 150 180];

% Gain points (in dB)
g1 = [6 3 2 1 .75 .5 .4 .3 .25 .2 .15 .1 .05 0 -.05 -.1 ...
      -.15 -.2 -.25 -.3 -.4 -.5 -.75 -1 -2 -3 -4 -5 -6 -9 ...
      -12 -16 -20:-10:max(-40,gmin) gmin(:,gmin<-40)];

% Compute gains GH and phases PH in H plane
[p,g] = meshgrid((pi/180)*p1,10.^(g1/20)); % in H/(1+H) plane
z = g .* exp(j*p);
H = z./(1-z);
gH = 20*log10(abs(H));
pH = rem((180/pi)*angle(H)+360,360);

% Add phase lines for angle between 180 and 360 (using symmetry)
gH = [gH , gH];
pH = [pH , 360-pH];

% Each column of GH/PH corresponds to one phase line
% Pad with NaN's and convert to vector
m = size(gH,2);
gH = [gH ; NanMat(1,ones(1,m))];
pH = [pH ; NanMat(1,ones(1,m))];
gain = gH(:);   phase = pH(:);


% (2) Generate isogain lines for following gain values:
g2 = [6 3 1 .5 .25 0 -1 -3 -6 -12 -20 -40:-20:gmin];

% Phase points
p2 = [1 2 3 4 5 7.5 10 15 20 25 30 45 60 75 90 105 ...
      120 135 150 175 180];
p2 = [p2 , fliplr(360-p2(1:end-1))];

[g,p] = meshgrid(10.^(g2/20),(pi/180)*p2);  % mesh in H/(1+H) plane
z = g .* exp(j*p);
H = z./(1-z);
gH = 20*log10(abs(H));
pH = rem((180/pi)*angle(H)+360,360);

% Each column of GH/PH describes one gain line
% Pad with NaN's and convert to vector
m = size(gH,2);
gH = [gH ; NanMat(1,ones(1,m))];
pH = [pH ; NanMat(1,ones(1,m))];
gain = [gain ; gH(:)];
phase = [phase ; pH(:)];

% Replicate Nichols chart if necessary
lp = length(phase);
dn = round((pmax-pmin)/360);     % number of 360 degree windup
gain = repmat(gain,dn,1);
phase = repmat(phase,dn,1);
shift = kron(pmin+360*(0:dn-1)',ones(lp,1));
ix = find(~isnan(phase));
phase(ix) = phase(ix) + shift(ix);
pH = pH + shift(end);

% Generate label output
labels = [pH(end-1,:) ; gH(end-1,:) ; g2];

