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
r9gmit.f
Go to the documentation of this file.
1 *DECK R9GMIT
2  FUNCTION r9gmit (A, X, ALGAP1, SGNGAM, ALX)
3 C***BEGIN PROLOGUE R9GMIT
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 SINGLE 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 ALNGAM, R1MACH, XERMSG
19 C***REVISION HISTORY (YYMMDD)
20 C 770701 DATE WRITTEN
21 C 890531 Changed all specific intrinsics to generic. (WRB)
22 C 890531 REVISION DATE from Version 3.2
23 C 891214 Prologue converted to Version 4.0 format. (BAB)
24 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
25 C 900720 Routine changed from user-callable to subsidiary. (WRB)
26 C***END PROLOGUE R9GMIT
27  SAVE eps, bot
28  DATA eps, bot / 2*0.0 /
29 C***FIRST EXECUTABLE STATEMENT R9GMIT
30  IF (eps.EQ.0.0) eps = 0.5*r1mach(3)
31  IF (bot.EQ.0.0) bot = log(r1mach(1))
32 C
33  IF (x .LE. 0.0) CALL xermsg('SLATEC', 'R9GMIT',
34  + 'X SHOULD BE GT 0', 1, 2)
35 C
36  ma = a + 0.5
37  IF (a.LT.0.0) ma = a - 0.5
38  aeps = a - ma
39 C
40  ae = a
41  IF (a.LT.(-0.5)) ae = aeps
42 C
43  t = 1.0
44  te = ae
45  s = t
46  DO 20 k=1,200
47  fk = k
48  te = -x*te/fk
49  t = te/(ae+fk)
50  s = s + t
51  IF (abs(t).LT.eps*abs(s)) go to 30
52  20 CONTINUE
53  CALL xermsg('SLATEC', 'R9GMIT',
54  + 'NO CONVERGENCE IN 200 TERMS OF TAYLOR-S SERIES', 2, 2)
55 C
56  30 IF (a.GE.(-0.5)) algs = -algap1 + log(s)
57  IF (a.GE.(-0.5)) go to 60
58 C
59  algs = -alngam(1.0+aeps) + log(s)
60  s = 1.0
61  m = -ma - 1
62  IF (m.EQ.0) go to 50
63  t = 1.0
64  DO 40 k=1,m
65  t = x*t/(aeps-m-1+k)
66  s = s + t
67  IF (abs(t).LT.eps*abs(s)) go to 50
68  40 CONTINUE
69 C
70  50 r9gmit = 0.0
71  algs = -ma*log(x) + algs
72  IF (s.EQ.0.0 .OR. aeps.EQ.0.0) go to 60
73 C
74  sgng2 = sgngam*sign(1.0,s)
75  alg2 = -x - algap1 + log(abs(s))
76 C
77  IF (alg2.GT.bot) r9gmit = sgng2*exp(alg2)
78  IF (algs.GT.bot) r9gmit = r9gmit + exp(algs)
79  RETURN
80 C
81  60 r9gmit = exp(algs)
82  RETURN
83 C
84  END