% programa 16.1, W.L.Luyben, pag. 566
% "multivariable Nyquist plot for Wood and Berry column"
% modif. E.Quintero.M 10nov'98
%-------------------------------------------------------
clear
format short e
% process transfer function data
kp=[12.8 -18.9; 6.6 -19.4];
n=2; D=[1 3; 7 3];
tau=[16.7 21; 10.9 14.4];
% empirical controller settings
%kc=[.2 -.04];
%reset=[4.44 2.67];
kc=[0.96 -0.19];
reset=[3.25 9.20];
%kc=[.2 -.04];
%reset=[0 0];
% set initial frequency and number of points per decade
w=0.01; dw=10^.1; xx=[];
% main loop for each frequency
while(w<10)
   s=i*w;
   for j=1:n
      for k=1:n
         gm(j,k)=kp(j,k)*exp(-D(j,k)*s)/(tau(j,k)*s+1);
      end
   end
   b=zeros(n);
   for j=1:n
      b(j,j)=kc(j)*(reset(j)*s+1)/(reset(j)*s); % quitar reset si igual a 0
   end
   q=gm*b;
   char=diag(ones(n,1),0)+q;
   xdeter=det(char);
   db=20*log10(abs(xdeter));
   wfunc=xdeter-1;
   deg=angle(wfunc)*180/pi;
   xx=[xx;w wfunc db deg];
w=w*dw;
end
plot(xx(:,2)),grid
axis([-2 2 -2 2])
