gammln.f

Go to the documentation of this file.

!                   ***************
                    FUNCTION gammln
!                   ***************
!
     &( xx    , deupi )
!
!***********************************************************************
! TOMAWAC   V6P1                                   15/06/2011
!***********************************************************************
!
!brief    COMPUTES THE NATURAL LOGARITHM FOR THE GAMMA FUNCTION
!+                (EULER FUNCTION OF SECOND-KIND).
!
!note     IF XX IS AN INTEGER NOTED N, GAMMA(N) = (N-1)!
!
!reference  "NUMERICAL RECIPES. THE ART OF SCIENTIFIC COMPUTING",
!+                       PRESS ET AL. (1989). (CF. PP 156-157)
!
!history  M. BENOIT
!+        15/11/95
!+        V1P0
!+   CREATED
!
!history  M. BENOIT
!+        07/11/96
!+        V1P2
!+   MODIFIED
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  G.MATTAROLO (EDF - LNHE)
!+        15/06/2011
!+        V6P1
!+   Translation of French names of the variables in argument
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| DEUPI          |-->| 2.PI
!| XX             |-->| VALUE AT WHICH LOG(GAMMA) IS CALCULATED
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_tomawac, ex_gammln => gammln
      IMPLICIT NONE
!
!.....VARIABLES IN ARGUMENT
!     """"""""""""""""""""
      DOUBLE PRECISION GAMMLN
      DOUBLE PRECISION,INTENT(IN)    ::  XX    , DEUPI
!
!.....LOCAL VARIABLES
!     """""""""""""""""
      INTEGER  J
      DOUBLE PRECISION STP   , X     , XC    , TMP   , SER   , AUX
      DOUBLE PRECISION COF(6)
!
!
      cof(1)= 76.180091730d0
      cof(2)=-86.505320330d0
      cof(3)= 24.014098220d0
      cof(4)= -1.231739516d0
      cof(5)=  0.001208580d0
      cof(6)= -0.000005364d0
      stp   =  2.506628275d0
!
      IF (xx.LT.1.d0) THEN
        xc=2.d0-xx
      ELSE
        xc=xx
      ENDIF
      x=xc-1.d0
      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp
      ser=1.d0
      DO j=1,6
        x=x+1.d0
        ser=ser+cof(j)/x
      ENDDO ! J
      gammln=tmp+log(stp*ser)
      IF (xx.LT.1d0) THEN
        aux=0.5d0*deupi*(1.d0-xx)
        gammln=log(aux/sin(aux))-gammln
      ENDIF
!
      RETURN
      END