C
C      **********************************************************************
C      *********************** FUNCTION PCHISQ  *****************************
C      **********************************************************************
C
       FUNCTION PCHISQ(CHISQR,NFREE)
C
C      *      THIS FUNCTION IS BASED ON PHILIP R. BEVINGTON 1969, "DATA     *
C      *      REDUCTION AND ERROR ANALYSIS FOR THE PHYSICAL SCIENCES", 1969,*
C      *      McGRAW HILL (NY:NY), PROGRAM 10-1 P. 192,                     *
C      *      AND COMPUTES CHI-SQUARE PROBABILITY FROM THE REDUCED          *
C      *      CHI-SQUARE.                                                   *
C      *      INPUT :      CHISQR : REDUCED CHI-SQUARE                      *
C      *                   NFREE  : DEGREES OF FREEDOM                      *
C      *     OUTPUT :      PCHISQ : CHI-SQUARE PROBABILITY                  *
C
       IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N)
C
       IF(NFREE .LE. 0) THEN
          PCHISQ=0.0
          RETURN
       ENDIF

       FREE=NFREE
       Z=CHISQR*FREE/2.0
       NEVEN=2*(NFREE/2)
       IF((NFREE-NEVEN) .LE. 0) THEN
C
C      *          THE DEGREES OF FREEDOM ARE EVEN                          *
C
   10     IMAX=NFREE/2
          TERM=1.0
          SUM=0.0
          DO 34 I=1,IMAX
             FI=I
             SUM=SUM+TERM
             TERM=TERM*Z/FI
   34     CONTINUE
          PCHISQ=SUM*DEXP(-Z)
          RETURN
C
C      *           THE DEGREES OF FREEDOM ARE ODD                          *
C
       ELSE
          IF((Z-25.0) .GT. 0) THEN
             Z=CHISQR*(FREE-1.0)/2.0
             IMAX=NFREE/2
             TERM=1.0
             SUM=0.0
             DO 44 I=1,IMAX
                FI=I
                SUM=SUM+TERM
                TERM=TERM*Z/FI
   44        CONTINUE
             PCHISQ=SUM*DEXP(-Z)
             RETURN
          ELSE
             PWR=FREE/2.0
             TERM=1.0
             SUM=TERM/PWR
             DO 56 I=1,1000
                FI=I
                TERM=-TERM*Z/FI
                SUM=SUM+TERM/(PWR+FI)
                IF((DABS(TERM/SUM)-0.00001) .LE. 0.0) GOTO 57
   56        CONTINUE
   57        PCHISQ=1.0-(Z**PWR)*SUM/GAMMA(PWR)
          ENDIF
       ENDIF

       RETURN
       END