# Help with creating my own Nyquist plotting function

7 views (last 30 days)
Daniel Beeson on 24 May 2020
Answered: Daniel Beeson on 24 May 2020
Hi everyone,
I am trying to make a function that takes a transfer function L(s) in the form of an array of numerator and denominator coefficients and spits out the nyquist plot, avoiding imaginary poles by drawing a semi-circle on the right half plane with radius epsilon around those imaginary poles and ending with an s-plane contour that is a semi-circle on the right half plane with radius R (see attached image)
So far, this is my code:
function n = nyquill(N,D,R,epsilon)
%N is numerator coefficients, D is denominator coefffs
L = tf(N,D);
t = pole(L)
r = NaN(3,1);
d = zeros(3,1);
for n = 1:size(t)
if real(t(n)) == 0
t(n) = r(n);
d(n) = NaN;
else
t(n) = d(n);
r(n) = NaN;
end
end
if r(1)==NaN & r(2)==NaN & r(3)==NaN
g1 = [-R:0.1:R]*i;
g2 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2];
elseif r(1)==NaN & r(2)==NaN
g1 = [-R;0.1:imag(r(3))-epsilon]*i;
g2 = imag(r(3)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(3)) + epsilon:.01:R]*i;
g4 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4];
elseif r(1)==NaN
g1 = [-R;0.1:imag(r(2))-epsilon]*i;
g2 = imag(r(2)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(2)) + epsilon:.01:imag(r(3))-epsilon]*i;
g4 = imag(r(3)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g5 = [imag(r(3)) + epsilon:.01:R]*i;
g6 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4 g5 g6];
elseif r(2)==NaN
g1 = [-R;0.1:imag(r(1))-epsilon]*i;
g2 = imag(r(1)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(1)) + epsilon:.01:imag(r(3))-epsilon]*i;
g4 = imag(r(3)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g5 = [imag(r(3)) + epsilon:.01:R]*i;
g6 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4 g5 g6];
elseif r(3)==NaN
g1 = [-R;0.1:imag(r(1))-epsilon]*i;
g2 = imag(r(1)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(1)) + epsilon:.01:imag(r(2))-epsilon]*i;
g4 = imag(r(2)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g5 = [imag(r(2)) + epsilon:.01:R]*i;
g6 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4 g5 g6];
elseif r(2)==NaN & r(3)==NaN
g1 = [-R;0.1:imag(r(1))-epsilon]*i;
g2 = imag(r(1)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(1)) + epsilon:.01:R]*i;
g4 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4];
elseif r(1)==NaN & r(3)==NaN
g1 = [-R;0.1:imag(r(2))-epsilon]*i;
g2 = imag(r(2)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(2)) + epsilon:.01:R]*i;
g4 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4];
elseif r(1)~=NaN & r(2)~=NaN & r(3)~=NaN
g1 = [-R;0.1:imag(r(1))-epsilon]*i;
g2 = imag(r(1)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(1)) + epsilon:.01:imag(r(2))-epsilon]*i;
g4 = imag(r(2)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g5 = [imag(r(2)) + epsilon:.01:imag(r(3))-epsilon]*i;
g6 = imag(r(3)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g7 = [imag(r(3)) + epsilon:.01:R]*i;
g8 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4 g5 g6 g7 g8];
end
c = poly2sym(N);
d = poly2sym(D);
ln = c/d;
x = g;
subs(ln);
plot(real(ln),imag(ln));grid on;axis('equal')
hold on
plot(-1,0, '-o')
end
And I end up receiving the error:
Error using plot
Data must be numeric, datetime, duration or an array convertible to double.
Error in nyquill (line 83)
plot(real(ln),imag(ln));grid on;axis('equal')
I don't know what this error means or exactly how to fix it, any help would be appreciated!

Daniel Beeson on 24 May 2020
The inputs N and D are the numerator and denominator coefficients
Daniel Beeson on 24 May 2020
The idea is to avoid the poles that lie on the imaginary axis by going around them with a small semi circle of radius epsilon that only goes into the right half plane

Daniel Beeson on 24 May 2020
Hey everyone, I figured it out. My conditions for NaN were incorrect. instead of using ==NaN, I changed it to isnan(r(...))==1 for if the value of r was NaN or ==0 for if it is not. Thank you everyone for the help!

rubindan on 24 May 2020