d9gmit.f

Go to the documentation of this file.

*DECK D9GMIT
      DOUBLE PRECISION FUNCTION d9gmit (A, X, ALGAP1, SGNGAM, ALX)
C***BEGIN PROLOGUE  D9GMIT
C***SUBSIDIARY
C***PURPOSE  Compute Tricomi's incomplete Gamma function for small
C            arguments.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7E
C***TYPE      DOUBLE PRECISION (R9GMIT-S, D9GMIT-D)
C***KEYWORDS  COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB, SMALL X,
C             SPECIAL FUNCTIONS, TRICOMI
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Compute Tricomi's incomplete gamma function for small X.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, DLNGAM, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   890911  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900720  Routine changed from user-callable to subsidiary.  (WRB)
C***END PROLOGUE  D9GMIT
      DOUBLE PRECISION a, x, algap1, sgngam, alx, ae, aeps, algs, alg2,
     1  bot, eps, fk, s, sgng2, t, te, d1mach, dlngam
      LOGICAL first
      SAVE eps, bot, first
      DATA first /.true./
C***FIRST EXECUTABLE STATEMENT  D9GMIT
      IF (first) THEN
         eps = 0.5d0*d1mach(3)
         bot = log(d1mach(1))
      ENDIF
      first = .false.
C
      IF (x .LE. 0.d0) CALL xermsg('SLATEC', 'D9GMIT',
     +   'X SHOULD BE GT 0', 1, 2)
C
      ma = a + 0.5d0
      IF (a.LT.0.d0) ma = a - 0.5d0
      aeps = a - ma
C
      ae = a
      IF (a.LT.(-0.5d0)) ae = aeps
C
      t = 1.d0
      te = ae
      s = t
      DO 20 k=1,200
        fk = k
        te = -x*te/fk
        t = te/(ae+fk)
        s = s + t
        IF (abs(t).LT.eps*abs(s)) go to 30
 20   CONTINUE
      CALL xermsg('SLATEC', 'D9GMIT',
     +   'NO CONVERGENCE IN 200 TERMS OF TAYLOR-S SERIES', 2, 2)
C
 30   IF (a.GE.(-0.5d0)) algs = -algap1 + log(s)
      IF (a.GE.(-0.5d0)) go to 60
C
      algs = -dlngam(1.d0+aeps) + log(s)
      s = 1.0d0
      m = -ma - 1
      IF (m.EQ.0) go to 50
      t = 1.0d0
      DO 40 k=1,m
        t = x*t/(aeps-(m+1-k))
        s = s + t
        IF (abs(t).LT.eps*abs(s)) go to 50
 40   CONTINUE
C
 50   d9gmit = 0.0d0
      algs = -ma*log(x) + algs
      IF (s.EQ.0.d0 .OR. aeps.EQ.0.d0) go to 60
C
      sgng2 = sgngam * sign(1.0d0, s)
      alg2 = -x - algap1 + log(abs(s))
C
      IF (alg2.GT.bot) d9gmit = sgng2 * exp(alg2)
      IF (algs.GT.bot) d9gmit = d9gmit + exp(algs)
      RETURN
C
 60   d9gmit = exp(algs)
      RETURN
C
      END