%  Example 5.3.  Nyquist plot generation.
clear
%  Define the plant transfer function.
nup=1;
dep=[1 -2 1];
%  Define the controller transfer function.
nuk=[1500 -100];
dek=[1 30 400];
%  Generate the loop transfer function.
nu=conv(nup,nuk);
de=conv(dep,dek);
%  Generate and plot the nyquist locus.
set(0,'DefaultAxesFontName','times')
set(0,'DefaultAxesFontSize',16)
set(0,'DefaultTextFontName','times')
set(0,'DefaultTextFontSize',16)
figure(1)
clf
w=logspace(-3,4,200);
[re,im]=nyquist(nu,de,w);
plot(re,im,'r-',re,-im,'r-',-1,0,'r+',re(75),im(75),...
   'r<',re(75),-im(75),'r>')
xlabel('Real axis')
ylabel('Imaginary axis')
axis([-2 2 -2 2])
axis('square')
grid
%  Add the unit circle onto the plot.
hold on
theta=0:.03:2*pi;
plot(cos(theta),sin(theta),'r--')
%  Add the line for measuring phase margin onto the plot.
x=[-1 0];
y=[-.98 0];
plot(x,y,'r--')
%  Add the phase angle arc onto the plot.
phi=-pi:.03:atan2(x(1),y(1));
plot(0.5*cos(phi),0.5*sin(phi),'r-')
%  Add the negative real axis onto the plot.
plot([-2 0],[0 0],'r-')
%  Add the labels: "-a", "-b", "PM" onto the plot
gtext('-a')
gtext('-b')
gtext('\it{\theta}')