GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
d9gmit.f
Go to the documentation of this file.
1 *DECK D9GMIT
2  DOUBLE PRECISION FUNCTION d9gmit (A, X, ALGAP1, SGNGAM, ALX)
3 C***BEGIN PROLOGUE D9GMIT
4 C***SUBSIDIARY
5 C***PURPOSE Compute Tricomi's incomplete Gamma function for small
6 C arguments.
7 C***LIBRARY SLATEC (FNLIB)
8 C***CATEGORY C7E
9 C***TYPE DOUBLE PRECISION (R9GMIT-S, D9GMIT-D)
10 C***KEYWORDS COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB, SMALL X,
11 C SPECIAL FUNCTIONS, TRICOMI
12 C***AUTHOR Fullerton, W., (LANL)
13 C***DESCRIPTION
14 C
15 C Compute Tricomi's incomplete gamma function for small X.
16 C
17 C***REFERENCES (NONE)
18 C***ROUTINES CALLED D1MACH, DLNGAM, XERMSG
19 C***REVISION HISTORY (YYMMDD)
20 C 770701 DATE WRITTEN
21 C 890531 Changed all specific intrinsics to generic. (WRB)
22 C 890911 Removed unnecessary intrinsics. (WRB)
23 C 890911 REVISION DATE from Version 3.2
24 C 891214 Prologue converted to Version 4.0 format. (BAB)
25 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
26 C 900720 Routine changed from user-callable to subsidiary. (WRB)
27 C***END PROLOGUE D9GMIT
28  DOUBLE PRECISION a, x, algap1, sgngam, alx, ae, aeps, algs, alg2,
29  1 bot, eps, fk, s, sgng2, t, te, d1mach, dlngam
30  LOGICAL first
31  SAVE eps, bot, first
32  DATA first /.true./
33 C***FIRST EXECUTABLE STATEMENT D9GMIT
34  IF (first) THEN
35  eps = 0.5d0*d1mach(3)
36  bot = log(d1mach(1))
37  ENDIF
38  first = .false.
39 C
40  IF (x .LE. 0.d0) CALL xermsg('SLATEC', 'D9GMIT',
41  + 'X SHOULD BE GT 0', 1, 2)
42 C
43  ma = a + 0.5d0
44  IF (a.LT.0.d0) ma = a - 0.5d0
45  aeps = a - ma
46 C
47  ae = a
48  IF (a.LT.(-0.5d0)) ae = aeps
49 C
50  t = 1.d0
51  te = ae
52  s = t
53  DO 20 k=1,200
54  fk = k
55  te = -x*te/fk
56  t = te/(ae+fk)
57  s = s + t
58  IF (abs(t).LT.eps*abs(s)) go to 30
59  20 CONTINUE
60  CALL xermsg('SLATEC', 'D9GMIT',
61  + 'NO CONVERGENCE IN 200 TERMS OF TAYLOR-S SERIES', 2, 2)
62 C
63  30 IF (a.GE.(-0.5d0)) algs = -algap1 + log(s)
64  IF (a.GE.(-0.5d0)) go to 60
65 C
66  algs = -dlngam(1.d0+aeps) + log(s)
67  s = 1.0d0
68  m = -ma - 1
69  IF (m.EQ.0) go to 50
70  t = 1.0d0
71  DO 40 k=1,m
72  t = x*t/(aeps-(m+1-k))
73  s = s + t
74  IF (abs(t).LT.eps*abs(s)) go to 50
75  40 CONTINUE
76 C
77  50 d9gmit = 0.0d0
78  algs = -ma*log(x) + algs
79  IF (s.EQ.0.d0 .OR. aeps.EQ.0.d0) go to 60
80 C
81  sgng2 = sgngam * sign(1.0d0, s)
82  alg2 = -x - algap1 + log(abs(s))
83 C
84  IF (alg2.GT.bot) d9gmit = sgng2 * exp(alg2)
85  IF (algs.GT.bot) d9gmit = d9gmit + exp(algs)
86  RETURN
87 C
88  60 d9gmit = exp(algs)
89  RETURN
90 C
91  END