function [re,im,s]=nyq(arg1,arg2,arg3,arg4,arg5,arg6,arg7)
%  [re,im,s]=nyq(num,den,w)
%  [re,im,s]=nyq(num,den,td,w)
%
%  Computes the Nyquist plot of open loop transfer function
%  L(s)= num(s)/den(s) using frequency vector w. The real and
%  imaginary parts of L(s) are returned separately for using
%  directly in plot(re,im). The vector  w need only be
%  the positive jw part of the Gamma contour normally used in 
%  the Nyquist plot, which 'encloses' the right half plane.
%  The small quarter-circle skipping the origin, the large 
%  quarter-circle enclosing the RHP and the mirror images of these
%  and the jw axis are automatically generated by NYQ, 
%  producing the complete Nyquist plot.  The optional td parameter
%  introduces a time lag in series with num(s)/den(s).
%  Example:
%
%  num=10;
%  den=[1 6 11 6];         i.e. L(s)=10/(s^3 +6s^2 +11s +6)
%  w=logspace(-2,4,100);
%  [re,im,s]=nyq(num,den,w)
%  plot(re,im)
%
%  [re,im,s]=nyq(A,B,C,D,iu,w)
%  [re,im,s]=nyq(A,B,C,D,td,iu,w)
%
%  computes the Nyquist plot of the state space realization
%  (A,B,C,D) from input number iu using frequency vector w.

%  Carlos Murillo, 2002/04/11
%  Took significant portions from the old bode.m function.
%  To do: compute margins w/o having to call bode.m .

nargs = nargin;
if ~((nargs==3)+(nargs==4)+(nargs==6)+(nargs==7))
  error('nyq: wrong number of input arguments');
end
if (nargs==3)
  num=arg1;
  den=arg2;
  w=arg3(:);
  td=0;
elseif (nargs==4)
  num=arg1;
  den=arg2;
  td=arg3;
  w=arg4;
elseif (nargs==6 | nargs==7)
  if nargs==6
    A=arg1;
    B=arg2;
    C=arg3;
    D=arg4;
    iu=arg5;
    w=arg6(:);
    td=0;
  else
    A=arg1;
    B=arg2;
    C=arg3;
    D=arg4;
    td=arg5;
    iu=arg6;
    w=arg7;
  end
  error(nargchk(6,7,nargs));
  error(abcdchk(A,B,C,D));
end

w=w(:);

% Generate small quarter circle z1 (avoiding open-loop poles in the origin)
% and large quarter-circle z2 (enclosing the farther upper-right-half-plane)
theta=(0:pi/16:pi/2)';
j=sqrt(-1);
z1=w(1)*exp(j*theta);
z2=w(length(w))*exp(j*theta);
z2=z2(length(z2):-1:1, :);

% Make one big vector with complex points enclosing the upper-RHP,
% i.e, half Nyquist contour.
z=[z1; j*w; z2];
% Now complete whole Nyquist contour
z = [ z ; conj(z(length(z):-1:1,:)) ];

% Eval. transfer function at complex contour
if (nargs == 3) | (nargs == 4)
	% It is in transfer function form.  Do directly, using Horner's method
	% of polynomial evaluation at the frequency points, for each row in
	% the numerator.  Then divide by the denominator.
        %[m,n] = size(num);
        %for i=1:m
        %g(:,i) = polyval(num(i,:),z);
        %end
        %g = (polyval(den,z)*ones(1,m).\g).';
        g = polyval(num,z) ./ polyval(den,z);
else
	[no,ns] = size(C);
	nz = max(size(z));

	% Balance A
	[t,A] = balance(A);
	B = t \ B;
	C = C * t;

	[p,A] = hess(A);		% Reduce A to Hessenberg form:

	%	 Apply similarity transformations from Hessenberg
	%	 reduction to B and C:
	B = p' * B(:,iu);
	C = C * p;
	D = D(:,iu);
	g = ltifr(A,B,z);
	g = C * g + diag(D) * ones(no,nz);
        g = g';
        g = g(length(g):-1:1);
end

% Introduce adjustments brought about by the time delay;
% convert g to column vector.
g = g .* exp(-z*td);
re_tmp = real(g);
im_tmp = imag(g);

if nargout == 0  % do plot only if no output arguments are detected
  ucirc = exp([ j*(0:0.02:2*pi),  0]);
  plot(re_tmp,im_tmp,real(ucirc),imag(ucirc), -1,0,'x', 0,0,'x', ...
     [-2 2],[0 0],'--',[0 0],[-2 2],'--')
  axis([-2 1 -1.5 1.5])
  axis square
else
  re = re_tmp;
  im = im_tmp;
  s = z;
end