% Displays the Popov plot for a dynamic linear system (integral action) H(s)
% with a time-invariant memoryless non-linear feedback controller phi(y) 
% residing in (0,k].
% Also outputs the largest k for which the closed-loop feedback
% system is absolutely stable.
% When prompted, enter the numerator and denominator coefficients for
% the transfer function H(s) without brackets.
% 
% Author : Abhishek Jain
% Date   : 10/11/2012

function popovplot
prompt = {'Enter numerator coefficients:',...
    'Enter denominator coefficients:'};
name = 'Popov plot';
numlines = 1;
defaultanswer = {'1','1 2 1 0'};
options.Resize='on';
options.WindowStyle='normal';
answer=inputdlg(prompt,name,numlines,defaultanswer,options);
[num,status] = str2num(answer{1});
if ~status
    num = 0;
end
[den,status] = str2num(answer{2});
if ~status
    den = 0;
end
h = tf(num,den);
w = logspace(-2,2,100000);
hw = freqs(num,den,w);
x = real(hw);
y = w.*imag(hw);
plot(x,y);
x0 = interp1(y,x,0);
max_k = -1/x0
grid on
xlabel('Re( H(jw) )')
ylabel('w * Im( H(jw) )')
T = evalc('h');
title(['Popov plot for',T])
text(x0,0.05,'-1/k')
hold on
plot(x0,0,'bo')
hold off
assignin('base','max_k',max_k);
end