v=gammac(z)

function v=gammac(z)
% v=gammac(z) gives the gamma function for complex argument.
% If no input is given, a surface plot is created.

%            HBW, August 17,2009
if nargin==0
  [x,y]=meshgrid(linspace(-5,4,80),linspace(-4,4,80));
  z=x+i*y;
end
c=1./[-122400/3617,156,-360360/691,1188,-1680,1260,-360,12];
siz=size(z); z=z(:); nz=length(z); logv=zeros(nz,1); N=10;
kb=find(real(z)>=N); zb=z(kb); nb=length(zb);
ks=find(real(z)<N); zs=z(ks); ns=length(zs); xmin=min(real(zs));
if xmin>0, n=N; else n=ceil(abs(xmin)+N); end
if nb>0
  logv(kb)=(zb-.5).*log(zb)-zb+.5*log(2*pi)+polyval(c,1./zb.^2)./zb;
end
if ns>0
  us=zs+n;  
  logv(ks)=(us-.5).*log(us)-us+.5*log(2*pi)+polyval(c,1./us.^2)./us;  
  logv(ks)=logv(ks)-sum(log(zs*ones(1,n)+ones(ns,1)*(0:n-1)),2);
end 
v=reshape(exp(logv),siz); 
if nargin==0
  close all  
  v=abs(v); v(find(v>=4))=4;  
  surf(x,y,abs(v)), colormap([1 1 0]),shg
  title('Absolute Value of the Gamma Function')
  xlabel('x axis'), ylabel('y axis'), zlabel('absolute(gamma(z))')
  v=[];
end