C GAMMA     SOURCE    KK2000    13/11/08    21:15:41     7860
      SUBROUTINE GAMMA(X,GA)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)

-INC CCREEL
C
C     ==================================================
C     Programme de : http://jin.ece.uiuc.edu/routines/routines.html
C     Purpose: Compute the gamma function G(x)
C     Input :  x  --- Argument of G(x)
C                     (x is not equal to 0,-1,-2,...)
C     Output:  GA --- G(x)
C     ==================================================
C
      DIMENSION G(26)
      DATA G/ 1D0                 , .5772156649015329D0,
     &       -.6558780715202538D0 ,-.420026350340952D-1,
     &        .1665386113822915D0 ,-.421977345555443D-1,
     &       -.96219715278770D-2  , .72189432466630D-2,
     &       -.11651675918591D-2  ,-.2152416741149D-3,
     &        .1280502823882D-3   ,-.201348547807D-4,
     &       -.12504934821D-5     , .11330272320D-5,
     &       -.2056338417D-6      , .61160950D-8,
     &        .50020075D-8        ,-.11812746D-8,
     &        .1043427D-9         , .77823D-11,
     &       -.36968D-11          , .51D-12,
     &       -.206D-13            ,-.54D-14,
     &        .14D-14             , .1D-15 /

C
      XINT=FLOAT(INT(X))
      IF (X.EQ.XINT) THEN
         IF (X.GT.0D0) THEN
            GA=1D0
            M1=INT(X-1.D0)
            DO 10 K=2,M1
10             GA=GA*FLOAT(K)
         ELSE
            GA=1D0+300D0
         ENDIF
      ELSE
         R=1D0
         IF (ABS(X).GT.1.D0) THEN
            Z=ABS(X)
            M=INT(Z)
            DO 15 K=1,M
15             R=R*(Z-FLOAT(K))
            Z=Z-FLOAT(M)
         ELSE
            Z=X
         ENDIF
         GR=G(26)
         DO 20 K=25,1,-1
20          GR=GR*Z+G(K)
         GA=1D0/(GR*Z)
         IF (ABS(X).GT.1.D0) THEN
            GA=GA*R
            IF (X.LT.0.D0) GA=-XPI/(X*GA*SIN(XPI*X))
         ENDIF
      ENDIF
      RETURN
      END