c-------------------------------------------------------------------
c     Bhattacharjee, G.P.(1970). 
c     The incomplete gamma integral 
c     Algorithm As 32, J.R.Statist.Soc.C, vol.19 no.3, pp. 285-287
c
c     available via the statlib archive, at http://lib.stat.cmu.edu/apstat/
c
c     function  gamain
c        computes incomplete gamma ratio for positive values of
c        arguments x and p.  g should equal ln(gamma(p)).
c        uses series expansion if p.gt.x or x.le.1, otherwise a
c        continued fraction approximation.
c
c     input variables:
c       x = upper limit of interal
c       p = power of integrand
c     output variables:
c       gamain = incomplete gamma ratio
c       ifault = 1 if p.le.0 else 2 if x.lt.0 else 0
c     auxiliary function required:
c       gamaln(p) = natural logarithm of p
c
c     last revision: 05-28-85 by jch
c-------------------------------------------------------------------
      doubleprecision function gamain(x,p,ifault)
      implicit double precision(a-h,o-z)
c6    implicit real*8 (a-h,o-z)
      dimension pn(6)
c
c       define accuracy and initialize
c
      acu=1.0d-8
      oflo=1.0d30
      gin=0.0d0
      ifault=0
c
c       test for admissibility of arguments
c
      if(p.le.0.0d0) ifault=1
      if(x.lt.0.0d0) ifault=2
      if(ifault.gt.0.or.x.eq.0.0d0) goto 50
      g=gamaln(p)
      factor=dexp(p*dlog(x)-x-g)
      if(x.gt.1.0d0.and.x.ge.p) goto 30
c
c       calculation by series expansion
c
      gin=1.0d0
      term=1.0d0
      rn=p
   20 rn=rn+1.0d0
      term=term*x/rn
      gin=gin+term
      if(term.gt.acu) goto 20
      gin=gin*factor/p
      goto 50
c
c       calculation by continued fraction
c
   30 a=1.0d0-p
      b=a+x+1.0d0
      term=0.0d0
      pn(1)=1.0d0
      pn(2)=x
      pn(3)=x+1.0d0
      pn(4)=x*b
      gin=pn(3)/pn(4)
   32 a=a+1.0d0
      b=b+2.0d0
      term=term+1.0d0
      an=a*term
      do 33 i=1,2
   33 pn(i+4)=b*pn(i+2)-an*pn(i)
      if(pn(6).eq.0.0d0) goto 35
      rn=pn(5)/pn(6)
      dif=dabs(gin-rn)
      if(dif.gt.acu) goto 34
      if(dif.le.acu*rn) goto 42
   34 gin=rn
   35 do 36 i=1,4
   36 pn(i)=pn(i+2)
      if(dabs(pn(5)).lt.oflo) goto 32
      do 41 i=1,4
   41 pn(i)=pn(i)/oflo
      goto 32
   42 gin=1.0d0-factor*gin
c
   50 gamain=gin
      return
      end
c