function makeNyquist(GH, w, ylimit, plotFilename)
fig = figure;
ylabel('Im(GH)')
grid on
xlim([-ylimit,ylimit])
ylim([-ylimit,ylimit])
axis equal
xlabel('Re(GH)')
hold on
line([0 0], ylim, 'color', 'k');
hold on
line(xlim, [0 0], 'color', 'k');
hold on
y = -1:0.001:1;
plot(sqrt(1-y.^2), y, '-k')
hold on
plot(-sqrt(1-y.^2), y, '-k')
hold on
plot(-1, 0, 'ok')
hold on
rgh = real(GH);
igh = imag(GH);
plot(rgh, igh, '-r')
hold on
plot(rgh, -igh, '--r')
hold on
for i = 2:length(igh)
if (igh(i)<0 && igh(i-1)>=0) || (igh(i)>0 && igh(i-1)<=0)
igh_xover_idx = i;
phz_xover_freq = w(igh_xover_idx);
reInt = rgh(igh_xover_idx);
GM = -20 * log10(abs(reInt));
sprintf('real axis intersect at Re = %.1f', reInt)
sprintf('gain margin (dB) = %.3f', GM)
plot(rgh(i), igh(i), '*b')
hold on
sprintf('phase crossover freq = %.3f', phz_xover_freq)
break
end
end
for j = 2:length(igh)
if (sqrt((rgh(j)).^2 + (igh(j)).^2)) <= 1
x = rgh(j);
y = igh(j);
plot(rgh(j), igh(j), '*b')
hold on
gain_xover_freq = w(j);
sprintf('gain crossover freq = %.3f', gain_xover_freq)
sprintf('unit circle intersect at (%.1f, %.1f)', x, y)
u = [x y];
v = [-1 0];
theta = acosd(dot(u,v) / (norm(u)*norm(v)));
if y > 0
theta = -theta;
elseif theta < 0
theta = theta + 180;
end
PM = theta;
sprintf('phase margin (deg) = %.3f', PM)
break
end
end
saveas(fig, plotFilename)
end