      FUNCTION pchi2(a,x)
      INTEGER i, iteven, ndeg
      REAL pchi2, a, x, gamma, prod, err, rest, f
      PARAMETER (err=5.e-7)
      EXTERNAL gamma

      IF (x.lt.0..or.a.le.0.) STOP 
     .     'Bad input for incomplete gamma function'
      ndeg=INT(a*2.)
      rest=MOD(a*2.,2.)
      
      IF (rest.lt.0.001) THEN ! even degrees of freedom
         iteven=ndeg/2
         pchi2=0.
         prod=1.
         DO i=1,iteven
            pchi2=pchi2+prod
            prod=prod*x/FLOAT(i)
         END DO
         pchi2=1.-EXP(-x)*pchi2
      ELSE                      ! odd degrees of freedom
         IF (x.ge.70.) THEN
            pchi2=1.
         ELSE
            c=1.
            pchi2=1.
            aa=a
 1          aa=aa+1.
            c=c*x/aa
            pchi2=pchi2+c
            IF (c/pchi2.gt.err) GO TO 1
            f=(x**a)*exp(-x)/gamma(a)
            pchi2=f*pchi2
         END IF
      END IF
      RETURN
      END

      FUNCTION qchi2(a,x)
      INTEGER i, iteven, ndeg
      REAL qchi2, a, x, gamma, prod, err, rest, f
      PARAMETER (err=5.e-7)
      EXTERNAL gamma

      IF (x.lt.0..or.a.le.0.) STOP 
     .     'Bad input for incomplete gamma function'
      ndeg=INT(a*2.)
      rest=MOD(a*2.,2.)
      
      IF (rest.lt.0.001) THEN ! even degrees of freedom
         iteven=ndeg/2
         qchi2=0.
         prod=1.
         DO i=1,iteven
            qchi2=qchi2+prod
            prod=prod*x/FLOAT(i)
         END DO
         qchi2=EXP(-x)*qchi2
      ELSE                      ! odd degrees of freedom
         IF (x.ge.80.) THEN
            qchi2=0.
         ELSE
            c=1.
            qchi2=1.
            aa=a
 1          aa=aa+1.
            c=c*x/aa
            qchi2=qchi2+c
            if (c/qchi2.gt.err) GO TO 1
            f=(x**a)*exp(-x)/gamma(a)
            qchi2=1.-f*qchi2
         END IF
      END IF
      RETURN
      END

      FUNCTION gamma(x)
      INTEGER i
      REAL x, gamma, xx
      REAL*8 coeff(6), fact
      DATA coeff/76.18009172947146D0,-86.50532032941677D0,
     &     24.01409824083091D0,-1.231739572450155D0,
     &     0.1208650973866179D-2,-0.5395239384953D-5/
      fact=1.000000000190015D0
      z=x
      DO i=1,6
         z=z+1.
         fact=fact+coeff(i)/z
      END DO
      xx=x+5.5
      gamma=xx**(x+0.5)*EXP(-xx)*SQRT(2.*ACOS(-1.))*fact
      RETURN
      END

      FUNCTION zero(func,xl,xu,err)
      REAL zero, func, xl, xu, x1, x2, xnew, err, test, errp
      INTEGER itmax
      PARAMETER (itmax=1000)
      EXTERNAL func
      IF (func(xl)*func(xu).gt.0.) THEN
         STOP ' ahi! zero non bracketed ...'
      END IF
      x1=xl
      x2=xu
      DO i=1,itmax
         xnew=x2-(func(x2)*(x1-x2))/(func(x1)-func(x2))
         IF (x1+x2.ne.0.) THEN
            errp=ABS((x2-x1)/(x2+x1))*100.
         END IF
         test=func(x1)*func(xnew)
         IF (test.eq.0.) THEN
            errp=0.
         ELSE
            IF (test.lt.0.) THEN
               x2=xnew
            ELSE
               x1=xnew
            END IF
         END IF
         IF (errp.le.err) GO TO 20
      END DO
 20   zero=xnew
      RETURN
      END
