function gl=lngamma(xx)
%LNGAMMA   Natural log of the complete Gamma function.
%	LNGAMMA(X) returns the log of the gamma function
%	for every element of X.  Accuracy is 1E-10.

% Peter R. Shaw, Woods Hole Oceanographic Institution
% Woods Hole, MA 02543
% (508) 457-2000 ext. 2473  pshaw@aqua.whoi.edu
% Revised GKS 5 June 92

% Converted from the Fortran subroutine "GAMMLN" in:
% Numerical Recipes, Press et al., Cambridge, 1986.

% ^ Calls no other routines ^

cof=[76.18009173, -86.50532033, 24.01409822, ...
-1.231739516, 0.120858003e-2, -0.536382e-5];
stp=2.50662827465;
[mxx,nxx]=size(xx);
one=ones(mxx,nxx);
half=0.5*one; fpf=5.5*one;
x=xx-one;
tmp=x+fpf;
tmp=(x+half).*log(tmp)-tmp;
ser=one;
for j=1:6,
  x=x+one;
  ser=ser+cof(j)./x;
end
gl=tmp+log(stp*ser);
