function computeNyquist(G, w)
%This function takes two inputs: Open loop transfer function G(s) as a
%symfun, and w as a vector of positive numbers

%let's first make a unit circle!! :D :D yay!
hold on

theta = 0:0.0001:2*pi
x_unit = cos(theta);
y_unit = imag(sin(theta).*1i); %y is the imaginary sin()
plot(x_unit,y_unit)

%making a little red circle at (-1 + 0j)
p1 = plot(-1,0,'or')
legend(p1, '-1 + j0')

%making the axis equal
axis equal
xlim([-2 2])
ylim([-2 2])
xlabel('Real Axis')
ylabel('Imaginary Axis')

%giving the plot a grid
grid on

%making some empty boxes for later
gain_marg = [];
phase_marg = [];
phase_xover_freq = [];
gain_xover_freq = [];

%plotting the positive frequencies

for a = 1:length(w)
   
    Gmag = double(abs(G(1i.*w(a))));
    p_angle = angle(G(1i.*w(a)));
    x = double(Gmag * cos(p_angle));
    y = double(imag(Gmag * sin(p_angle) * 1i));
    if a == length(w)
        p2 = plot(x,y,'*b') %the last point is a blue star and represents +infinity
        legend(p2, '+infinity')
    else
        plot(x,y,'.g') %all other points are green dots
    end

    if a > 1
        
        Gmag_prev = double(abs(G(1i.*w(a-1))));
        p_angle_prev = angle(G(1i.*w(a-1)));
        y_prev = double(imag(Gmag_prev * sin(p_angle_prev) * 1i));
        
        if (y_prev >= 0 && y <= 0) || (y_prev <= 0 && y >= 0) %if plot crosses the real axis 
            phase_xover_freq = w(a);
            gain_marg = -20*log10(abs(x));    
        end
        if (Gmag_prev >= 1 && Gmag <= 1) || (Gmag_prev <= 1 && Gmag >= 1)
            gain_xover_freq = w(a);
            if p_angle <= 0
                phase_marg = double(p_angle * 180/pi + 180);
            else
                phase_marg = double(p_angle * 180/pi - 180);
            end
        end

    end

end

w = -1 .* w; % now for the negative frequencies

for a = 1:length(w)
   
    Gmag = abs(G(1i.*w(a)));
    p_angle = angle(G(1i.*w(a)));
    x = Gmag * cos(p_angle);
    y = imag(Gmag * sin(p_angle) * 1i);
    if a == length(w)
        p3 = plot(x,y,'*m') %the last point is a magenta star and represents -infinity
    else
        plot(x,y,'.y') %all other points are yellow dots
    end
end

legend([p1 p2 p3], '-1 + j0', 'towards +infinity', 'towards -infinity')

clc

if ~isempty(gain_marg)
    sprintf('The gain margin is %0.4g dB.', round(gain_marg,2))
else
    sprintf('The gain margin is undefined.')
end

if ~isempty(phase_marg)
    sprintf('The phase margin is %0.4g degrees.', round(phase_marg,2))
else
    sprintf('The phase margin is undefined.')
end

if ~isempty(gain_xover_freq)
    sprintf('The gain crossover freq is %0.4g radians per second.', round(gain_xover_freq,2))
else
    sprintf('The gain crossover freq is undefined.')
end

if ~isempty(phase_xover_freq)
    sprintf('The phase crossover freq is %0.4g radians per second.', round(phase_xover_freq,2))
else
    sprintf('The phase crossover freq is undefined.')
end

end
