function zpplot(zepo,sd,mode,ax)
%ZPPLOT  Plots zeros and poles.
%
%	zpplot(ZEPO)   or   zpplot(ZEPO,SD)  or  zpplot(ZEPO,MODE)
%		       or   zpplot(ZEPO,SD,MODE,AXIS)
%
%	The zeros and poles, specified by ZEPO (See TH2ZP, or ZPFORM 
%	for the format) are plotted, with 'o' denoting zeros and 'x' poles.
%	Poles and zeros	associated with the same input, but different models
%	are always plotted in the same diagram, and 'ENTER' advances the plot 
%	from one model to the next (if any).
%
%	When ZEPO contains information about several different inputs there
%	are some options: 
%	MODE='sub' (The default value) splits the screen into several plots.
%	MODE='same' gives all plots in the same diagram. Use 'ENTER' to advance
%	MODE='sep' erases the previous plot before the next input is treated.
%       AXIS=[x1 x2 y1 y2] fixes the axes scaling.
%	AXIS=m is the same as AXIS=[-m m -m m]. Default is autoscaling
%	If SD>0, confidence regions around the poles and zeros are plotted.
%	The region corresponding to SD standard deviations is marked.
%	SD=0 is default.
%	For multi-output models, each output is handled separately.

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

% Some intial tests on the input arguments
if nargin<4, ax=[];end
if nargin<3, mode=[];end
if nargin<2,sd=0;end, if isempty(sd),sd=0;end,if sd<0, sd=0;end
if length(sd)>1 mode=sd; sd=0;end
if isempty(ax),axdef=1; else axdef=0;end
if isempty(mode),mode='sub';end
mode=mode(1:3);
if mode=='sub',sbflag=1;else sbflag=0;end
[nzp,nc]=size(zepo);
if fix(nc/2)~=nc/2
  error('ZEPO should have an even number of columns')
end

%
colors=['g','b','r','w'];
zsc=zepo;
if any(any(zsc==inf));ztemp=zsc(:);
ztemp(find(ztemp==inf))=zeros(length(find(ztemp==inf)),1);zsc(:)=ztemp;end
om=2.1*pi*[1:100]/100;
w=exp(om*sqrt(-1));

inputs1=zepo(1,:);
if abs(sum(sign(inputs1)))~=nc
disp(['WARNING: ZEPO contains a mixture of discrete and continuous models.'...
'\nPress a key to continue, C-c to abort']),pause,end
T=sign(inputs1(1));inputs1=abs(inputs1);zepo(1,:)=inputs1;
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

for ko=kko+1
inputs=inputs1-(ko-1)*1000;

kki=[]; mi=max(inputs(find(inputs<20)));% The input indices
mino=max(inputs(find(inputs>500 & inputs<520)));
for ki=[1:mi 501:mino]
if length(find(inputs==ki))>0,kki=[kki ki];end
end

clf reset
space=[];
sbcount = 0;
if mode=='sub' & norm(kki-kki(1)*ones(1,length(kki)))>0,
   sbcount=1;end
if mode=='sam',mm=max(max(abs(zsc(2:nzp,find((inputs)/60<1)))));
    m=ceil(mm);if T<0,mar=max(max(real(zsc(2:nzp,find((inputs)/60<1)))));
    mir=min(min(real(zsc(2:nzp,find((inputs)/60<1)))));end
    
end
for ki=kki
    inpc=find((inputs==ki) | (inputs==ki+20));
    insd=find((inputs==ki+60) | (inputs==ki+80));
    if sum(mode~='sam')>0,
       mm=max(max(abs([zsc(2:nzp,inpc)]))); m=ceil(mm);
       if T<0,mar=max(max(real([zsc(2:nzp,inpc)]))); 
       mir=min(min(real([zsc(2:nzp,inpc)]))); end
    end


    if ~axdef m=ax; end
    if length(m)==1, m1=m;m(1)=-m1;m(2)=m1;m(3)=-m1;m(4)=m1;end
    if sbflag & sbcount>0,subplot(220+sbcount),end
    plot([m(1);m(2)],[m(3);m(4)],'i.'),axis(axis)
    axis('square'),

    
    hold on
    if sum(mode~='sam')>0
       if ki<500 ,title(['OUTPUT # ',int2str(ko),'  INPUT # ',int2str(ki)])
       else title(['OUTPUT # ',int2str(ko),'  NOISE INPUT # ',int2str(rem(ki,500))])
       end
    else text(0.3,1.03,['OUTPUT # ',int2str(ko),' INPUT # ',space,int2str(ki)],'sc')
    space=[space,'  '];
    end

    if T>0,plot(w),end %The unit circle
    count=0;
    for k=inpc
    	colind=rem(count/2,4)+1;count=count+1;
        if k+1>nc,sdpl=0;info=zepo(1,k)-1000*(ko-1);
	if info>500,info=info-500;end
	else info=zepo(1,k:k+1)-1000*(ko-1);if info(1,1)>500,info=info-500;end
	sdpl=((info(1,2)==info(1,1)+60)&sd>0);
	end
	if sdpl
	   zp=zepo(2:nzp,k:k+1);
	else zp=zepo(2:nzp,k);
	end
	iz=find(abs(zp(:,1))<1/eps);
	if T>0,iz=find((abs(zp(:,1))>0).*(abs(zp(:,1))<1/eps));end
	if length(iz)>0
          if fix(info(1,1)/20)==0,mark1='o';mark2='-.';mark3='-';
	  else mark1='x';mark2=':';mark3='-';
	  end
	  mark1=[mark1,colors(colind)];
	  mark2=[mark2,colors(colind)];
	  mark3=[mark3,colors(colind)];
          plot(real(zp(iz,1)),imag(zp(iz,1)),mark1)
	  if sdpl,zpsdpl(zp,sd,w,iz',mark2,mark3),end
        end
	if fix(info(1,1)/20)>0 & k~=inpc(length(inpc)), pause,end
    end
    if sbflag,sbcount=sbcount+1;
    if sbcount>4,sbcount=1;end,end
    if sum(mode~='sam')>0 ,hold off,end
    if ki~=kki(length(kki)),pause,end
end
if ko~=kko(length(kko))+1,pause,hold off,end
end

hold off

set(gcf,'NextPlot','replace');
