function show_zser(a, opt);
% function show_zser(a, opt);
%
% a = [a0 a1 a2 ... an] power series coeffs (if empty plots opt.func)
% opt is options struct, all its fields are optional:
% opt.ngrid : # gridpoints in x and y directions
% opt.win   : viewing window half-width on complex plane
% opt.z0    : expansion point (if not 0)
% opt.func  : function to compare against, eg '1./(1-z)' (note ./)
% opt.slice : if nonzero, relative height of real slice plot, eg 0.3)
%
% Alex Barnett 11/11/05, some based on code of Matthias Kawski.

if ~isfield(opt,'ngrid')
  opt.ngrid = 100;
end
if ~isfield(opt,'win')
  opt.win = 5;
end
if ~isfield(opt,'z0')
  opt.z0 = 0 + 0i;
end
grid = -opt.win : 2*opt.win/opt.ngrid : opt.win;
[x,y] = meshgrid(grid, grid);
z = x + i*y;
if isempty(a)
  f = eval(opt.func);
else
  f = polyval(fliplr(a), z-opt.z0);
end
if isfield(opt,'slice')
  clf
  axes('position', [.05 opt.slice+.05 .9 0.9-opt.slice]);
end
p = surf(x, y, zeros(size(x)), 'CData', z2rgb_kawski(f));
xlabel('Re[z]'); ylabel('Im[z]');
set(p,'LineStyle','none','FaceColor','interp');
view(2);
axis equal
if isfield(opt,'slice')
axes('position', [.07 .07 .84 opt.slice-.07]);
  z = real(opt.z0)+grid;
  plot(z, real(polyval(fliplr(a), z-opt.z0)), '-');
  xlabel('Re[z]'); ylabel('Re[f]');
  if isfield(opt,'func')
    fx = real(eval(opt.func));
    hold on; plot(z, fx, 'g--'); legend(sprintf('n=%d',numel(a)-1), opt.func);
    axis([min(z) max(z) max(-5,min(fx)) min(5,max(fx))]);
  end
end

function [rgb] = z2rgb_kawski(z)
% function [rgb] = z2rgb_kawski(z)
%
% converts complex numbers into values in unit RGB cube
% 
% Alex Barnett 11/11/05, influenced by Matthias Kawski
m = abs(z) + 1e-15;
s = 1 - (m/2 + 1).^-1;
x = real(z)./m;           % x,y is point on unit circle
y = imag(z)./m;
w = 3*s.*(1-s);   % radius of cone as func of saturation, adjust prefactor
rgb(:,:,1) = s + w.*(-y/sqrt(6) + x/sqrt(2));
rgb(:,:,2) = s + w.*(-y/sqrt(6) - x/sqrt(2));
rgb(:,:,3) = s + w.*y*sqrt(2/3);