      function gamma_p(a,x) ! berg june 10 1999.  
c incomplete gamma function. 
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      parameter (iter_max=3800,eps=3.d-9)
      lbad1=.false.
      lbad2=.false.
      if(x.lt.zero.or.a.le.zero) print*,'a,x:',a,x
      if(x.lt.zero.or.a.le.zero) stop 'gamma_p - input not covered'
      gln=gamma_ln(a)
3     continue
c     print*,"3 continue:"
      if(x.lt.a+one.and..not.lbad2.or.lbad1) then ! series expansion
c       print*,"1. a,x:",a,x
        gamma_p=zero
        if(x.le.zero) return
        sum=one/a
        add=sum
        do iter=1,iter_max
          add=add*x/(a+iter)
          sum=sum+add
        if(abs(add).lt.abs(sum)*eps) goto 1
        end do 
        if(lbad1) stop 'gamma_p-1: a too large, iter_max too small'
        lbad2=.true.
      go to 3
1       gamma_p=sum*exp(-x+a*log(x)-gln)
c       print*,"1. gamma_p:",gamma_p
      else ! continued fraction expansion
c       print*,"2. a,x,gamma_p:",a,x,gamma_p
        g=zero
        a0=x+one-a
        b0=one
        a1=x*a0+x
        b1=x*b0+one
        do iter=2,iter_max
          gold=g
          xiter=iter*one
          xiterma=xiter-a
          a0=(a1+a0*xiterma)/a1
          b0=(b1+b0*xiterma)/a1
          b1=x*b0+xiter*b1/a1
          a1=x*a0+xiter
          if(a1.eq.zero) stop 'a1 zero in gamma_p'
          g=b1/a1
          if(abs((g-gold)/g).lt.eps) go to 2
        end do
      if(lbad2) stop 'gamma_p-2: a too large, iter_max too small'
        lbad1=.true.
      go to 3
2       gamma_p=one-exp(-x+a*log(x)-gln)*g ! and take its complement
c       print*,"2. iter,gln,g,gamma_p:",iter,gln,g,gamma_p
      endif
      return
      end
