Help with creating my own Nyquist plotting function

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]

The inputs N and D are the numerator and denominator coefficients

[Daniel Beeson]

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

Answer 1

[Daniel Beeson]

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!

Answer 2

[rubindan]

See nyquist1.m