incog Subroutine

subroutine incog ( a, x, gin, gim, gip )

  !*****************************************************************************80
  !
  !! INCOG computes the incomplete gamma function r(a,x), Γ(a,x), P(a,x).
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    22 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, real ( kind = 8 ) A, the parameter.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) GIN, GIM, GIP, the values of
  !    r(a,x), Γ(a,x), P(a,x).
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) ga
  real ( kind = 8 ) gim
  real ( kind = 8 ) gin
  real ( kind = 8 ) gip
  integer ( kind = 4 ) k
  real ( kind = 8 ) r
  real ( kind = 8 ) s
  real ( kind = 8 ) t0
  real ( kind = 8 ) x
  real ( kind = 8 ) xam

  xam = -  x + a * log ( x )

  if ( 700.0D+00 < xam .or. 170.0D+00 < a ) then
     write ( *, '(a)' ) ' '
     write ( *, '(a)' ) 'INCOG - Fatal error!'
     write ( *, '(a)' ) '  A and/or X is too large!'
     stop
  end if

  if ( x == 0.0D+00 ) then

     gin = 0.0D+00
     call gammaf ( a, ga )
     gim = ga
     gip = 0.0D+00

  else if ( x <= 1.0D+00 + a ) then

     s = 1.0D+00 / a
     r = s
     do k = 1, 60
        r = r * x / ( a + k )
        s = s + r
        if ( abs ( r / s ) < 1.0D-15 ) then
           exit
        end if
     end do

     gin = exp ( xam ) * s
     call gammaf ( a, ga )
     gip = gin / ga
     gim = ga - gin

  else if ( 1.0D+00 + a < x ) then

     t0 = 0.0D+00
     do k = 60, 1, -1
        t0 = ( k - a ) / ( 1.0D+00 + k / ( x + t0 ) )
     end do
     gim = exp ( xam ) / ( x + t0 )
     call gammaf ( a, ga )
     gin = ga - gim
     gip = 1.0D+00 - gim / ga

  end if

  return
end subroutine incog