bief_gcf.f

Go to the documentation of this file.

!                       *******************
                        SUBROUTINE bief_gcf
!                       *******************
!
     &(gammcf,a,x,gln)
!
!***********************************************************************
! BIEF   V6P3                                   07/03/2013
!***********************************************************************
!
!brief  Returns the incomplete gamma function Q(a,x) evaluated by its
!+      continued fraction representation as GAMMCF.  Also returns
!+      gamma (a) as GLN. From 'Numerical Recipes' Chapter 6.2
!
!history  D J Evans-Roberts (HRW)
!+        23/11/1993
!+        V6P3
!+   First version sent by Michiel Knaapen (HRW) on 07/03/2013.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A              |-->| PARAMETER
!| GAMMCF         |<->| Incomplete gamma function Q(a,x)
!| GLN            |<->| gamma(a)
!| X              |-->| OPERAND
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     USE BIEF
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN)    :: A,X
      DOUBLE PRECISION, INTENT(INOUT) :: GAMMCF,GLN
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER :: N
      DOUBLE PRECISION GOLD,A0,A1,B0,B1,FAC,AN,ANA,ANF,G
!
      INTEGER, PARAMETER :: ITMAX = 100
      DOUBLE PRECISION   :: EPS = 3.e-7
!
      DOUBLE PRECISION BIEF_GAMMLN
      EXTERNAL         bief_gammln
!
      INTRINSIC exp,log,float,abs
!
!-----------------------------------------------------------------------
!
      gln = bief_gammln(a)
!     This is the previous value, tested against for convergence.
      gold = 0.d0
      a0   = 1.d0
!     We are here setting up the A's and B's of equation (5.2.4) for
!     evaluating the continued fraction.
      a1     = x
      b0     = 0.d0
      b1     = 1.d0
!     FAC is the renormalization factor for preventing overflow of the
!     partial numerators and denominators.
      fac    = 1.d0
      DO n = 1,itmax
        an     = float(n)
        ana    = an - a
!       One step of the recurrence (5.2.5).
        a0     = (a1+a0*ana)*fac
        b0     = (b1+b0*ana)*fac
        anf    = an*fac
!       The next step of the recurrence (5.2.5).
        a1     = x*a0 + anf*a1
        b1     = x*b0 + anf*b1
!       Shall we renormalize?
        IF(a1.NE.0.d0) THEN
!         Yes.  Set FAC so it happens.
          fac    = 1.d0/a1
!         New value of answer.
          g      = b1*fac
!         Converged? If so, exit.
          IF (abs((g-gold)/g).LT.eps) GO TO 110
!         If not, save value.
          gold   = g
        ENDIF
      ENDDO
      WRITE(lu,*) 'BIEF_GCF: FATAL ERROR'
      CALL plante(1)
      stop
110   CONTINUE
!
      gammcf = exp(-x+a*log(x)-gln)*g
!
!-----------------------------------------------------------------------
!
      RETURN
      END