bief_gser.f

Go to the documentation of this file.

!                       ********************
                        SUBROUTINE bief_gser
!                       ********************
!
     &(gamser,a,x,gln)
!
!***********************************************************************
! BIEF   V6P3                                   07/03/2013
!***********************************************************************
!
!brief  Returns the incomplete gamma function P(a,x) evaluated by its
!+      series representation as GAMSER.  Also returns ln  gamma(a) as
!+      GLN. From 'Numerical Recipes' Chapter 6.2 This routine is only
!+      valid for X.GE.0 Negative X cannot occur when calling from ERF
!+      or ERFC via GAMMP or GAMMQ
!
!history  D J Evans-Roberts (HRW)
!+        23/11/1993
!+        V6P3
!+   First version sent by Michiel Knaapen (HRW) on 07/03/2013.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A              |-->| PARAMETER
!| GAMSER         |<->| Incomplete gamma function P(a,x)
!| GLN            |<->| gamma(a)
!| X              |-->| OPERAND
!| IFAIL          |<->| ERROR TAG: 1 IS OK, 0 IS ERROR
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     USE BIEF
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN)    :: X,A
      DOUBLE PRECISION, INTENT(INOUT) :: GAMSER,GLN
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER :: N
      DOUBLE PRECISION :: AP,SOMME,DEL
!
      INTEGER, PARAMETER :: ITMAX = 100
      DOUBLE PRECISION   :: EPS = 3.d-7
!
      DOUBLE PRECISION BIEF_GAMMLN
      EXTERNAL         bief_gammln
!
      INTRINSIC abs,exp,log
!
!-----------------------------------------------------------------------
!
      gln = bief_gammln(a)
!
      IF(x.LT.0.d0) THEN
!       Cannot occur when used by ERF or ERFC
        WRITE(lu,*) 'BIEF_GSER: UNEXPECTED NEGATIVE X'
        CALL plante(1)
        stop
      ELSEIF(x.EQ.0.d0) THEN
        gamser = 0.d0
      ELSE
        ap  = a
        somme = 1.d0/a
        del = somme
        DO n = 1,itmax
          ap  = ap + 1.d0
          del = del*x/ap
          somme = somme + del
          IF(abs(del).LT.abs(somme)*eps) GO TO 110
        ENDDO
        WRITE(lu,*) 'BIEF_GSER: MAXIMUM OF ITERATIONS REACHED:',itmax
        CALL plante(1)
        stop
110     CONTINUE
        gamser = somme*exp(-x+a*log(x)-gln)
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END