function show_z2(B,A);
%function show_z2(B,A);
%
% SHOW_Z2 --  a function that displays a surface plot of the Z-plane
%             for a given filter defined by B and A.  Also displayed
%             are the magnitude frequency response and a pole-zero plot
%
%  Of course, the z-plane "blows up" in many locations, so this function
%  "cuts off" the surface at 2+ the maximum point on the unit circle.


num_pts = 46;
upplus = 2;

[B1,A1] = eqtflength(B,A);

[z,p,k] = tf2zp(B1,A1);

mn = min([-1.2; real([z; p]); imag([z;p])]);
mx = max([1.2; real([z; p]); imag([z;p])]);

[x,y] = meshgrid(linspace(mn,mx,num_pts),linspace(mn,mx,num_pts));

z = x+i*y;

%zi = 1./z;

ce = exp(-j*2*pi*(0:.001:1));
cevals = abs(polyval(B,ce)./polyval(A,ce));

vals = polyval(B(end:-1:1),1./z)./polyval(A(end:-1:1),1./z);

%vals(abs(vals) > max(cevals)+upplus) = NaN;
vals(abs(vals) > max(cevals)+upplus) = max(cevals)+upplus;

subplot(121);
mesh(x,y,abs(vals));

axis tight;

h = line(real(ce),imag(ce),cevals + 0.01);
set(h,'Color',[1 0 0]);
set(h,'LineWidth',2);

h = line(real(ce),imag(ce),cevals - 0.01);
set(h,'Color',[1 0 0]);
set(h,'LineWidth',2);

xlabel('Real(z)');
ylabel('Imag(z)');
zlabel('|H(z)|');

subplot(222);
zplane(B,A);

subplot(224);
[H,w] = freqz(B,A);
plot(w,abs(H),'r');

xlabel('Discrete radian frequency, \omega');
ylabel('|H(\omega)|');
axis([min(w),max(w),0,max(abs(H))]);
grid on
