c$Header: /weru/cvs/weps/weps.src/util/math/gammln.for,v 1.1.1.1 1999-03-12 17:05:31 wagner Exp $
c
      real function gammln
     i                 (xx)
c
c     + + + PURPOSE + + +
c     Computes the ln of the gamma function for xx > 0
c     Full accuracy is obtained for xx > 1
c     For 0 < xx < 1, the reflection formula (6.1.4) can be used first.
c
c     Based on:
c     'NUMERICAL RECIPES - The Art of Scientific Computing',
c     W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling
c     Cambridge University Press, 1986
c     pg 157
c
c     + + + KEYWORDS + + +
c     gamma function
c
c     + + + ARGUMENT DECLARATIONS + + +
      real   xx
c
c
c     + + + ARGUMENT DEFINITIONS + + +
c     xx     - real value for values > 0
c
c     + + + LOCAL VARIABLES + + +
      double precision cof(6), stp
      double precision half, one, fpf
      double precision x, tmp, ser
      integer j
c
      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
c
c     + + + END SPECIFICATIONS + + +
c
      x = dble(xx)-one
      tmp = x+fpf
      tmp = (x+half)*log(tmp)-tmp
      ser=one
      do 100 j=1,6
          x = x+one
          ser = ser+cof(j)/x
100   continue
      gammln = tmp+log(stp*ser)

      return
      end