Draw the Nichols grid

% Chapter 15 Exercise 14.
% — Draw the Nichols grid
xmin = -360; xmax = 0; ymin=-20; ymax=20;
margins =[xmin xmax ymin ymax];
ha=axis(margins);
degs =xmin:5:xmax;
dbs  = ymin:0.5:ymax;
[DEG, DB] = meshgrid(degs, dbs);
ol    = 10.^(DB/20).*                …
exp(sqrt(-1)*DEG*pi/180);      % open loop, numerical value
cl    = ol./(1+ol);                    % closed loop
cldb  = 20*log10(abs(cl));             % closed loop, decibel
vals  = [-12 -6 0 6 12];               % values to plot
% — plot the M-circles
contour(DEG,DB,cldb,vals);
% — plot tics and rectangular grid
xticks=xmin:45:xmax;                   % x-axis tics
yticks=ymin:5:ymax;                    % y-axis ticks
set(gca,’XTick’,xticks);
set(gca,’YTick’,yticks);
grid on;
hold on;

% — Draw the transfer function of the plant
pG     = 5.0;                          % plant gain
pP     = 10;                           % plant pole
pzi    = 0.7;                          % plant damping factor
pomn   = 15;                           % plant natural frequency
pnum   = [0 0 0 0 pG*pomn^2];
pden   = conv([1 0], conv([1/pP 1],  …
[1 2*pzi*pomn pomn^2]));

om     = logspace(-1, 2);
s      = sqrt(-1)*om;

P      = polyval(pnum,s) ./ polyval(pden,s);
magdb  = 20*log10(abs(P));
phadeg = 180/pi*unwrap(angle(P));
plot(phadeg, magdb);

disp(‘ To determine gain margin, read the point where the graph    ‘);
disp(‘ intersects the vertical line of -180 degrees.               ‘);
disp(‘ The gain margin is then the opposite than the value you read’);
disp(‘ for the magnitude                                           ‘);
disp(”);
disp(‘ Press any key when ready …                                ‘);
pause
[gpha, gdb]=ginput(1);
disp([‘Gain margin: ‘, num2str(-gdb), ‘ db’]);

disp(‘ To determine phase margin, read the point where the graph   ‘);
disp(‘ intersects the horizontal line of 0 db.                     ‘);
disp(‘ The phase margin is then 180 plus the value you read for the’);
disp(‘ phase.                                                      ‘);
disp(”);
disp(‘ Press any key when ready …                                ‘);
pause

[pdeg, pdb]=ginput(1);
disp([‘Phase margin: ‘, num2str(180+pdeg), ‘ deg’]);

disp(‘ If you have access to the function MARGIN from the Control  ‘);
disp(‘ Toolbox, you can verify your answer against that computed   ‘);
disp(‘ by MATLAB, otherwise, you will get an error message.        ‘);
disp(‘ Press any key when ready …                                ‘);
pause
[Gm, Pm, Wcg, Wcp]=margin(pnum, pden);
disp([‘MATLAB reports      Gain margin:  ‘,num2str(20*log10(Gm))]);
disp([‘                    Phase margin: ‘,num2str(Pm)])