function zplane3(zeros,poles);
%ZPLANE	Z-plane zero-pole plot using 3-D coordinates
%	ZPLANE3(Z,P) mesh plot of the pole-zero pattern including a
%   unit circle for reference.  
%
%   If both inputs are column vectors, they are assumed to be
%   zeros and poles.
%	If either argument is a row, ZPLANE3 finds the roots of the argument
%	with ROOTS.  So, ZPLANE3(B,A) where B and A are row vectors containing
%	transfer function polynomial coefficients plots the poles and zeros
%	of B(z)/A(z).  ZPLANE3 assumes scalars are zeros or poles.
%
%	See also FREQZ.

%	Author(s): T. Krauss, 3-19-93
%	Copyright (c) 1984-94 by The MathWorks, Inc.
%       $Revision: 1.11 $  $Date: 1994/01/25 18:00:11 $

%   Modified by Richard Stern, 3/11/97, 9/30/97, 9/20/99

	newplot;
    if (nargin<2),
        p=[];
    end;
   
    if (prod(size(zeros))==size(zeros,2))     % convert row vector to zeros
        if prod(size(zeros))>1   %  (but not if it's a scalar)
            zeros = roots(zeros);
        end
    end
    if (prod(size(poles))==size(poles,2))     % convert row vector to poles
        if prod(size(poles))>1   %  (but not if it's a scalar)
            poles = roots(poles);
        end
    end

    ax = gca;

    if ~any(imag(zeros)),
        zeros = zeros + j*1e-50;
    end;
    if ~any(imag(poles)),
        poles = poles + j*1e-50;
    end;
 
realz = -1.5:.05:1.5;
imagz = -1.5:.05:1.5;
[Realz Imagz] = meshgrid(realz,imagz);
z = Realz+j*Imagz;

% Now set up the unit circle:

w = 0:pi/50:2*pi;
realw = cos(w);
imagw = sin(w);
zw = realw +j* imagw;


% Now calculate the response

if length(zeros) > 0
	b = poly(zeros);
else
	b = 1;
end

if length(poles) > 0
	a = poly(poles);
else
	a = 1;
end

amax = length(a);
bmax = length(b);
imax = max(amax,bmax);

Hnum = b(bmax);
Hwnum = b(bmax);
Hdenom = a(amax);
Hwdenom = a(amax);

ztemp = z;
zwtemp = zw;

for i = 2:imax,
	if(i <= bmax),
		Hnum = Hnum + (b(bmax-i+1).*ztemp);
		Hwnum = Hwnum + (b(bmax-i+1).*zwtemp);
		end
	if(i <= amax),
		Hdenom = Hdenom + (a(amax-i+1).*ztemp);
		Hwdenom = Hwdenom + (a(amax-i+1).*zwtemp);
		end
		
ztemp = z.*ztemp;
zwtemp = zw.*zwtemp;
end

H = Hnum./(eps+Hdenom);
Hw = Hwnum./(eps+Hwdenom);

% mesh(realz,imagz,abs(H))		
mesh(realz,imagz,min(10*max(abs(Hw)),abs(H)))		% Clamp max to +60 dB
hidden off
hold on
plot3(realw,imagw,abs(Hw),'r')
%axis([-1.5 1.5 -1.5 1.5 0 10])
hold off
view(30,45)
xlabel('Real[z]')
ylabel('Imag[z]')