function bodeplot(G,sd,C,mode)
%BODEPLOT  plots the Bode diagram of a transfer function or spectrum
%
%	bodeplot(G)  or  bodeplot(G,SD)
%
%	G is entered as the transfer function(s) or spectrum (spectra)
%	of the standard frequency function format, produced by TH2FF,
%	SPA or ETFE. See HELP FREQFUNC.
%	
%	Several plots in the same diagram are obtained by
%	bodeplot([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 amplitude and phase plots are shown simul-
%	taneously, for each input (noise spectrum) present in G. For spectra
%	the phase plots are omitted. ENTER will advance the plot from one
%	input to the next (if any). To show amplitude or phase plots only 
%	use bodeplot(G,SD,'A') and bodeplot(G,SD,'P'), respectively. To ob-
%	tain all plots in the same diagrams use bodeplot(G,SD,C,'same').
%	C is here 'A', 'P' or 'B' for amplitude, phase or both plots.
%	(In these cases SD can be omitted for default)

%	L. Ljung 10-1-86, revised 7-7-87,10-2-92
%	Copyright (c) 1986-92 by the MathWorks, Inc.
%	All Rights Reserved.

if nargin<4 , mode=[];end
if nargin<3, C=[];end
if nargin<2, sd=[];end
if isempty(mode),mode='sep';end
if isempty(C),C='B';end
if isempty(sd),sd=0;end

abp=['a';'A';'b';'B';'p';'P'];
if sum(sd(1)==abp)>0, if length(C)>1,mode=C;end,C=sd;sd=0;end
if length(sd)>1,mode=sd;sd=0;end

mode=mode(1:3);
if mode=='SEP', mode='sep';end
if mode=='SAM', mode='sam';end
if mode~='sam', mode='sep';end

if C == 'a', C='A'; end
if C == 'p', C='P'; end
if C == 'b', C='B'; end

[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 sum(kki)==0 & C=='B', C='A';end
if mode=='sep'
  for k=kki
  if k==0 & C=='B',C='A';end
  clf reset,if C=='B',subplot(211),else subplot(111),end
  if C=='A'| C=='B'

	loglog(G(:,find(inputs==100+k)),G(:,find(inputs==k)))
	if sd>0,hold on,
	ind=find(k+50==inputs);
	loglog(G(:,ind-2),G(:,ind-1)+sd*G(:,ind),'-.')
	loglog(G(:,ind-2),G(:,ind-1)-sd*G(:,ind),'-.') %Change here 2.04
	hold off,end
	xlabel('frequency (rad/sec)')

	if k~=0,title(['AMPLITUDE PLOT, input # ',int2str(k),' output # ',int2str(ko)])
	else title(['SPECTRUM output # ',int2str(ko)]),end
  end
  if ((C=='B')|(C=='P'))&k~=0
        if (C=='B'), subplot(212), end
        if (C=='P'),clf,subplot(111),end

	semilogx(G(:,find(inputs==100+k)),G(:,find(inputs==20+k)))
	if sd>0,hold on,
	ind=find(k+70==inputs);
	semilogx(G(:,ind-4),G(:,ind-1)+G(:,ind)*sd,'-.')
	semilogx(G(:,ind-4),G(:,ind-1)-sd*G(:,ind),'-.') %Change here 2.04
	hold off,end
	title(['PHASE PLOT, input # ',int2str(k),' output # ',int2str(ko)])
	xlabel('frequency (rad/sec)')
	ylabel('phase')
  end


  if k~=kki(length(kki)),pause,end
end
end
if mode=='sam' 

	if C=='B',subplot(211),else subplot(111),end
      if C=='B' | C=='A'	
	loglog(G(:,find(inputs>99&inputs<120)),G(:,find(inputs<20&inputs>-1)))
	if sd>0,hold on,
	ind=find(inputs<69 & inputs>49);
	loglog(G(:,ind-2),G(:,ind-1)+G(:,ind)*sd,'-.')
	loglog(G(:,ind-2),G(:,ind-1)-sd*G(:,ind),'-.'),hold off,end
	xlabel('frequency (rad/sec)')
	title(['AMPLITUDE PLOT output #',int2str(ko)])
      end
      if C=='B'| C=='P'
         if C=='B',subplot(212),else subplot(111),end
        semilogx(G(:,find(inputs>100&inputs<120)),G(:,find(inputs<39 & inputs>19)))
	if sd>0,hold on
	ind=find(inputs<89 & inputs>69);
	semilogx(G(:,ind-4),G(:,ind-1)+G(:,ind)*sd,'-.')
	semilogx(G(:,ind-4),G(:,ind-1)-sd*G(:,ind),'-.'),hold off,end
	xlabel('frequency (rad/sec)')
	ylabel('phase')
	title(['PHASE PLOT output #',int2str(ko)])
      end
	
end,if ko~=kko(length(kko))+1,pause,end

end
hold off
set(gca,'NextPlot','replace');
set(gcf,'NextPlot','replace');
