function nyqplot(G,sd,mode)
%NYQPLOT  plots the Nyquist diagram of a transfer function 
%
%	nyqplot(G)  or  nyqplot(G,SD)
%
%	G is entered as the transfer function(s) of the standard frequency 
%	function format, produced by TH2FF, SPA or ETFE. 
%	See HELP FREQFUNC.
%	
%	Several plots in the same diagram are obtained by
%	nyqplot([G1 G2 .. Gn]) 
%	The transfer functions Gk need not be given at the same frequencies.
%	If SD is specified to a number larger than zero, also confidence
%	intervals for the functions (corresponding to SD standard deviations)
%	are shown as dash-dotted  curves.
%
%	The default mode is that ENTER will advance the plot from one
%	input to the next (if any). To obtain all plots in the same diagram
%	use nyqplot(G,SD,'same').
%	Frequency functions for different outputs will always be given in
%	separate plots.

%	L. Ljung 10-7-89
%	Copyright (c) 1989-90 by the MathWorks, Inc.
%	All Rights Reserved.

if (nargin==1), sd=0; mode='sep';end
if (nargin==2),mode='sep';end
clf
mode=mode(1:3);

[m,n]=size(G);
inputs1=G(1,:);mo=fix(inputs1/1000);mmo=max(mo);
kko=[]; % The output indices
for k=[0:mmo],if length(find(mo==k))>0,kko=[kko k];end,end
G=G(2:m,:);
for ko=kko+1
inputs=inputs1-(ko-1)*1000;
mi=max(inputs(find(inputs<19)));

kki=[]; % The input indices
for k=[1:mi 0], if length(find(inputs==k))>0,kki=[kki k];end,end

if mode=='sep'
for k=kki
if k~=0,      
clg

	gp=G(:,find(inputs==20+k))*pi/180;
	polar(gp,G(:,find(inputs==k)))
	ax=axis; hold on, polar([0 0],ax(1:2)),
	polar([pi pi]/2,ax(1:2)),hold off
	if sd>0,hold on,
	ind=find(k+50==inputs);
	polar(gp,G(:,ind-1)+sd*G(:,ind),'-.')
	polar(gp, G(:,ind-1)-sd*G(:,ind),'-.')
	hold off,end

	

	title(['NYQUIST PLOT, input # ',int2str(k),' output # ',int2str(ko)])
	
else disp('No Nyquist plot for spectra!'),end


if k~=kki(length(kki)),pause,end
end
end

if mode=='sam' 
		gp=G(:,find(inputs<39 & inputs > 19))*pi/180;
      	polar(gp,G(:,find(inputs<20)))
	ax=axis; hold on, polar([0 0],ax(1:2)),
	polar([pi pi]/2,ax(1:2)),hold off
	if sd>0,hold on,
        ind=find(inputs<69 & inputs>49);
	polar(gp,G(:,ind-1)+sd*G(:,ind),'-.')
	polar(gp, G(:,ind-1)-sd*G(:,ind),'-.')
	
	hold off,end
	
	title(['NYQUIST PLOT output # ',int2str(ko)])
      end
if ko~=kko(length(kko))+1,pause,end
end
hold off
set(gcf,'NextPlot','replace');

