%
% Macro to produce complex root-locus plots
% (c) Marc Bodson, June 18, 2020
%
% Inputs:
%    num, den: row vector polynomials with deg(num) <= deg(den)
%    gmax: max gain for the root-locus
%    npt: number of points per branch
%
% Pre-computations
%
g=(0:1:npt-1)*gmax/(npt-1);
numext=[zeros(1,length(den)-length(num)) num];
nb=length(den)-1;ar=zeros(npt,nb);
%
% Compute the root-locus
%
for i=1:npt
    poly=den+g(i)*numext;
    ar(i,:)=roots(poly);
end
%
% Connect the branches
%
for i=2:npt
    ind=1:nb;v1=ar(i-1,:)';v2=ar(i,:)';
    v2out=zeros(size(v1));
    for ib=1:nb
        dv=abs(v1(ib)-v2(ind));
        [dvmin,dvi]=min(dv);
        v2out(ib)=v2(ind(dvi));
        if ib<nb;ind(dvi)=[];end
    end
    ar(i,:)=v2out';
end
%
% Plot the root-locus
%
cols=['b' 'g' 'r' 'c' 'm' 'y' 'k'];ncols=length(cols);
figure(1);clf(1);
hold on;
for ib=1:nb
    plot(real(ar(:,ib)),imag(ar(:,ib)),cols(mod(ib-1,ncols)+1));
end
%
% Mark the open-loop poles and zeros
%
pol=roots(den);zol=roots(num);
for i=1:length(pol)
    plot(real(pol(i)),imag(pol(i)),'xk','MarkerSize',12)
end
for i=1:length(zol)
    plot(real(zol(i)),imag(zol(i)),'ok','MarkerSize',10)
end
xlabel('Re(s)')
ylabel('Im(s)')
hold off