bief_gammln.f

Go to the documentation of this file.

!                 *************************************
                  DOUBLE PRECISION FUNCTION bief_gammln
!                 *************************************
!
     &(xx)
!
!***********************************************************************
! BIEF   V6P3                                   07/03/2013
!***********************************************************************
!
!brief    Returns the value ln(gamma(xx) for xx > 0) Full accuracy is
!+        obtained for xx > 1.
!
!history  D J Evans-Roberts (HRW)
!+        23/11/1993
!+        V6P3
!+   First version sent by Michiel Knaapen (HRW) on 07/03/2013.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| XX             |-->| OPERAND
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     USE BIEF
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN) :: XX
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER J
      DOUBLE PRECISION :: COF(6),STP,HALF,ONE,FPF,X,TMP,SER
      parameter( cof = (/ 76.18009173d0,-86.50532033d0,24.01409822d0,
     &             -1.231739516d0,0.120858003d-2,-0.536382d-5 /) )
      parameter( stp = 2.506628727465d0 )
      parameter( half = 0.5d0, one = 1.d0, fpf = 5.5d0 )
!
      INTRINSIC log
!
!-----------------------------------------------------------------------
!
      x   = xx - one
      tmp = x + fpf
      tmp = (x+half)*log(tmp) - tmp
      ser = one
!
      DO j = 1,6
        x   = x + one
        ser = ser + cof(j)/x
      ENDDO
!
      bief_gammln = tmp + log(stp*ser)
!
!-----------------------------------------------------------------------
!
      RETURN
      END