Complex plane of LTI model visualization

Below is the octave/matlab function cplane which will visualize complex transfer function of an LTI model with a single input and a single output (SISO) into a 3D mesh or surface. The LTI model can be either a laplace 's' domain or discrete 'z' domain model. The file is available here.
% USAGE:
% cplane (SYS, RES, AX, FUNC, RESP)
% Plot the magnitude or phase of a transfer function in a complex plane
% Inputs:
% 	SYS
% 		LTI model, must be a single input and single output model
% 	RES 
% 		(optional) Real and imaginary resolution of complex plane
% 	AX
% 		(optional) Axis limits of real and imaginary axes.
% 	FUNC
% 		(optional) define mesh plot function. Options are: 'mesh', 'meshc','meshz', 'surf', 'surfc', 'pcolor'
%	RESP
%		(optional) specify type of response, either "mag" (default) or "phase".
% Output:
%	[oxx,oyy,ozz]
%		(optional) The coordinatates to call plot functions, or export to other tools
%
% See also pzmap, bode, bodemag
%
% examples:
%	% 4th order butterworth low-pass filter, corner frequency=1000 [rad/s]
%	[num, den]=butter(4,1000,'s'); 
%	h=tf(num,den);
%	cplane(h);
%
%	%lowpass elliptic/Cauer filter, 5th-order,  0.1/20dB's pass/stop ripple, fc=20 [rad/s]
%	[num, den]=ellip(5,0.1,20,20,'s'); 
%	h=tf(num,den);
%	cplane(h,150,[-20 0 -30 30]);
%
%	%bessel 9th order low-pass filter, fc=100 [rad/s]
%	[num, den]=besself(9,100,'s');
%	h=tf(num,den);
%	cplane(h,[],[], 'surfc'); % use the contour surface function
%	view(0,90); shading interp; colormap gray; %create a log density plot view
%

% Copyright (C) 2012 Ronald van der Meer
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2,
% or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; If not, see <http://www.gnu.org/licenses/>.
function [oxx, oyy, ozz] = cplane(sys,res,ax,func,resp)
  %narginchk(1, 5) %check the correct number of arguments
  if (nargin < 5) || (exist('resp','var')==0) || (length(resp)==0);
    resp='mag'; %default response type
  end
  if (nargin < 4) || (exist('func','var')==0) || (length(func)==0);
    func='mesh'; %default plot function
  end
  if (nargin < 3) || (exist('ax','var')==0) || (length(ax)==0); %determine if auto-scaling is needed
    pz=[pole(sys) ;zero(sys)];
    rmin=min(real(pz));
    rmax=max(real(pz));
    imin=min(imag(pz));
    imax=max(imag(pz));
    a=0.2; %This is the scaling factor e.g. 0.2 = 20% beyond the min/max of the pole/zero's
    ax=[rmin-(rmax-rmin)*a rmax+(rmax-rmin)*a imin-(imax-imin)*a imax+(imax-imin)*a];
  end
  if (nargin < 2) || (exist('res','var')==0) || (length(res)==0);
    res=80; %the default xy resolution of the plot
  end;
  x=linspace(ax(1),ax(2),res);
  y=linspace(ax(3),ax(4),res);
  [xx,yy]=meshgrid(x,y);
  z=polyval(sys.num{1},xx+i*yy)./polyval(sys.den{1},xx+i*yy);
  if (strcmp(resp,'mag'))
    zz=20*log10(abs(z)); %convert to [dB] magnitude
  else
    zz=atan2(imag(z),real(z)); %do phase instead of magnitide
  end
  if (nargout==0) %when there are no output arguments
    eval(strcat(func,'(x,y,zz)'));
    axis('tight');
    set(gca,'TickDir','out') %let the tick-lines be outside the plot. This will position the tick-text a little better
    xlabel('Real Axis');
    ylabel('Imaginary Axis');
    if (strcmp(resp,'mag'))
      zlabel('Magnitude [dB]');
    else
      zlabel('Phase [rad]');
    end
    %change the view so we look "into" the negetive real half plane
    [az,el]=view;
    az=az+90;
    view(az,el);
  else %don't plot but return coordinates
    oxx=x;
    oyy=y;
    ozz=zz;
  end
end

%!demo
%! [num, den]=ellip(5,0.1,20,20,'s'); 
%! h=tf(num,den);
%! cplane(h);

Here is some example octave code that produces a 4th order Butterworth filter complex plane plot:

graphics_toolkit fltk
figure
[num, den]=butter(4,1000,'s'); 
h=tf(num,den);
cplane(h);
print -dpng butter4example.png -FArialMT:8 "-S800,600"
system('convert butter4example.png -trim butter4example.png')

Here is some example octave code that produces an elliptic 5-th order filter complex plane plot:

graphics_toolkit fltk
figure
[num, den]=ellip(5,0.1,20,10,'s');  
h=tf(num,den);
cplane(h,150,[],'meshc');
print -dpng cauer5example.png -FArialMT:8 "-S1200,1000"
system('convert cauer5example.png -trim cauer5example.png')

Here is some example octave code for a pole-zero map, with matching log density plot:

graphics_toolkit fltk
figure
[num, den]=ellip(5,0.1,20,10,'s'); 
h=tf(num,den);
subplot(1,2,1)
pzmap(h);
legend('poles','zeros');
axislimits=[get(gca,'xlim') , get(gca,'ylim')]
axis('equal');
subplot(1,2,2);
cplane(h,80,axislimits);
set(gca,'TickDir','in') %let the tick-lines be inside to match pzmap
title('Log density of h')
axis('equal');
view(0,90); shading interp; colormap gray;
print -dpng logdensityexample.png -FArialMT:8 "-S1000,1100"
system('convert logdensityexample.png -trim logdensityexample.png')

Here is some example octave code for a notch filter with matching frequency magnitude and phase plots:

%graphics_toolkit fltk
figure(1)
wc=10; %center freq.
bw=10; %-3dB bandwidth=10, Q-factor=1
num = [1 0 wc^2];  %notch numerator polynomials
den = [1 bw wc^2]; %notch denominator polynomials
h=tf(num,den);
cplane(h,[],[-6 0 0 16 ]);
print -dpng notchexample_tmp1.png -FArialMT:6 "-S540,350"
system('convert notchexample_tmp1.png -trim notchexample_cmag.png')

figure(2)
cplane(h,[],[-6 0 0 16 ],[],'phase')
print -dpng notchexample_tmp2.png -FArialMT:6 "-S540,350"
system('convert notchexample_tmp2.png -trim notchexample_cphase.png')

w=linspace(0,16,100);
[mag,phase,w]=bode(h,w);
figure(3)
plot(w,20*log(mag))
xlabel('freq in [rad/s]');
ylabel('Magnitude [dB]');
grid;
axis tight;
print -dpng notchexample_tmp3.png -FArialMT:6 "-S490,200"
system('convert notchexample_tmp3.png -trim notchexample_mag.png')

figure(4)
plot(w,phase/180*pi)
xlabel('freq in [rad/s]');
ylabel('Phase [rad]');
grid;
axis tight;
print -dpng notchexample_tmp4.png -FArialMT:6 "-S490,200"
system('convert notchexample_tmp4.png -trim notchexample_phase.png')