C
C      **********************************************************************
C      ********************* FUNCTION GAMMA  ********************************
C      **********************************************************************
C
       FUNCTION GAMMA(X)
C
C      *     THIS FUNCTION IS OBTAINED FROM PHILIP R. BEVINGTON, "DATA      *
C      *     REDUCTION AND ERROR ANALYSIS FOR THE PHYSICAL SCIENCES", 1969, *
C      *     McGRAW HILL (NY:NY), PROGRAM 7-2 P. 126                        *
C      *     THIS COMPUTES THE GAMMA FUNCTION FOR INTEGERS AND HALF-INTEGERS*
C
       IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N)
C
C      *                INTEGERIZE ARGUMENT                                 *
C
       N = INT(X - 0.25)
       XN=REAL(N)
   13  IF((X-XN-0.75) .GT. 0.0) THEN 
C
C      *                ARGUMENT IS INTEGER                                 *
C
          GAMMA=FACTOR(N)
C
C      *                ARGUMENT IS HALF-INTEGER                            *
C
       ELSE
   31     PROD=1.77245385D00
          IF(N .LE. 0) THEN
             GAMMA=PROD
             GOTO 56
          ENDIF
          IF(N-10 .LE. 0) THEN
             DO 43 I=1,N
                FI=REAL(I)
                PROD=PROD*(FI-0.5)
   43        CONTINUE
   44        GAMMA=PROD
          ELSE
   51        SUM=0.0
             DO 54 I=11,N
                FI=I
                SUM=SUM+DLOG(FI-0.5)
   54        CONTINUE
   55        GAMMA=PROD*639383.8623D00*DEXP(SUM)
          ENDIF
   56  ENDIF
   60  RETURN
       END