r9gmit.f

Go to the documentation of this file.

*DECK R9GMIT
      FUNCTION r9gmit (A, X, ALGAP1, SGNGAM, ALX)
C***BEGIN PROLOGUE  R9GMIT
C***SUBSIDIARY
C***PURPOSE  Compute Tricomi's incomplete Gamma function for small
C            arguments.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7E
C***TYPE      SINGLE 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  ALNGAM, R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  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  R9GMIT
      SAVE eps, bot
      DATA eps, bot / 2*0.0 /
C***FIRST EXECUTABLE STATEMENT  R9GMIT
      IF (eps.EQ.0.0) eps = 0.5*r1mach(3)
      IF (bot.EQ.0.0) bot = log(r1mach(1))
C
      IF (x .LE. 0.0) CALL xermsg ('SLATEC', 'R9GMIT',
     +   'X SHOULD BE GT 0', 1, 2)
C
      ma = a + 0.5
      IF (a.LT.0.0) ma = a - 0.5
      aeps = a - ma
C
      ae = a
      IF (a.LT.(-0.5)) ae = aeps
C
      t = 1.0
      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', 'R9GMIT',
     +   'NO CONVERGENCE IN 200 TERMS OF TAYLOR-S SERIES', 2, 2)
C
 30   IF (a.GE.(-0.5)) algs = -algap1 + log(s)
      IF (a.GE.(-0.5)) GO TO 60
C
      algs = -alngam(1.0+aeps) + log(s)
      s = 1.0
      m = -ma - 1
      IF (m.EQ.0) GO TO 50
      t = 1.0
      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   r9gmit = 0.0
      algs = -ma*log(x) + algs
      IF (s.EQ.0.0 .OR. aeps.EQ.0.0) GO TO 60
C
      sgng2 = sgngam*sign(1.0,s)
      alg2 = -x - algap1 + log(abs(s))
C
      IF (alg2.GT.bot) r9gmit = sgng2*exp(alg2)
      IF (algs.GT.bot) r9gmit = r9gmit + exp(algs)
      RETURN
C
 60   r9gmit = exp(algs)
      RETURN
C
      END