c  Andrew M. O'Reilly 11/15/2001
c
c*****************************************************
c  Gamma Function, approximation from:
c  Press, W.H., Flannery, B.P., Teukolsky, S.A., and
c       Vetterling, W.T., 1989, Numerical recipes --
c       The art of scientific computing (FORTRAN version):
c       Cambridge University Press, 702 p.
c
c  The error of this numerical approximation is reported
c  to be less than 2e-10 for xx>=1, and for 0<xx<1 when
c  the suggested reflection formula is used.
c
c  This version is for xx>0 and uses the reflection
c  formula for 0<xx<1 for additional accuracy.
c
c  The natural log of the gamma function is used to
c  avoid overflow.
c******************************************************
c
      function gammln(xx)
      implicit double precision (a-h,o-z)
	dimension cof(6)
	data cof,stp /76.18009173d0,-86.50532033d0,24.01409822d0,
     &-1.231739516d0,0.120858003d-2,-0.536382d-5,2.50662827465d0/
	data half,one,fpf /0.5d0,1.0d0,5.5d0/
c
	if (xx.gt.0.and.xx.lt.1.) then	!use reflection formula
	   irefl=1
	   pi=4.d0*datan(1.d0)
	   x=one-xx
	elseif (xx.ge.1.) then   !do not use reflection formula
	   irefl=0
	   x=xx-one
	else
	   write(*,*) 'INVALID--Gamma Function argument must be > 0'
	endif
c
	tmp=x+fpf
	tmp=(x+half)*log(tmp)-tmp
	ser=one
	do 10 j=1,6
	   x=x+one
	   ser=ser+cof(j)/x
   10 continue
      gammln=tmp+log(stp*ser)
c
	if (irefl.eq.1) then
	   gammln=log(pi*(one-xx))-gammln-log(sin(pi*(one-xx)))
	endif
c
	return
	end