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
gamit.f
Go to the documentation of this file.
1 *DECK GAMIT
2  REAL FUNCTION gamit (A, X)
3 C***BEGIN PROLOGUE GAMIT
4 C***PURPOSE Calculate Tricomi's form of the incomplete Gamma function.
5 C***LIBRARY SLATEC (FNLIB)
6 C***CATEGORY C7E
7 C***TYPE SINGLE PRECISION (GAMIT-S, DGAMIT-D)
8 C***KEYWORDS COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB,
9 C SPECIAL FUNCTIONS, TRICOMI
10 C***AUTHOR Fullerton, W., (LANL)
11 C***DESCRIPTION
12 C
13 C Evaluate Tricomi's incomplete gamma function defined by
14 C
15 C GAMIT = X**(-A)/GAMMA(A) * integral from 0 to X of EXP(-T) *
16 C T**(A-1.)
17 C
18 C for A .GT. 0.0 and by analytic continuation for A .LE. 0.0.
19 C GAMMA(X) is the complete gamma function of X.
20 C
21 C GAMIT is evaluated for arbitrary real values of A and for non-
22 C negative values of X (even though GAMIT is defined for X .LT.
23 C 0.0), except that for X = 0 and A .LE. 0.0, GAMIT is infinite,
24 C which is a fatal error.
25 C
26 C The function and both arguments are REAL.
27 C
28 C A slight deterioration of 2 or 3 digits accuracy will occur when
29 C GAMIT is very large or very small in absolute value, because log-
30 C arithmic variables are used. Also, if the parameter A is very
31 C close to a negative integer (but not a negative integer), there is
32 C a loss of accuracy, which is reported if the result is less than
33 C half machine precision.
34 C
35 C***REFERENCES W. Gautschi, A computational procedure for incomplete
36 C gamma functions, ACM Transactions on Mathematical
37 C Software 5, 4 (December 1979), pp. 466-481.
38 C W. Gautschi, Incomplete gamma functions, Algorithm 542,
39 C ACM Transactions on Mathematical Software 5, 4
40 C (December 1979), pp. 482-489.
41 C***ROUTINES CALLED ALGAMS, ALNGAM, GAMR, R1MACH, R9GMIT, R9LGIC,
42 C R9LGIT, XERCLR, XERMSG
43 C***REVISION HISTORY (YYMMDD)
44 C 770701 DATE WRITTEN
45 C 890531 Changed all specific intrinsics to generic. (WRB)
46 C 890531 REVISION DATE from Version 3.2
47 C 891214 Prologue converted to Version 4.0 format. (BAB)
48 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
49 C 920528 DESCRIPTION and REFERENCES sections revised. (WRB)
50 C***END PROLOGUE GAMIT
51  LOGICAL first
52  SAVE alneps, sqeps, bot, first
53  DATA first /.true./
54 C***FIRST EXECUTABLE STATEMENT GAMIT
55  IF (first) THEN
56  alneps = -log(r1mach(3))
57  sqeps = sqrt(r1mach(4))
58  bot = log(r1mach(1))
59  ENDIF
60  first = .false.
61 C
62  IF (x .LT. 0.0) CALL xermsg('SLATEC', 'GAMIT', 'X IS NEGATIVE',
63  + 2, 2)
64 C
65  IF (x.NE.0.0) alx = log(x)
66  sga = 1.0
67  IF (a.NE.0.0) sga = sign(1.0, a)
68  ainta = aint(a+0.5*sga)
69  aeps = a - ainta
70 C
71  IF (x.GT.0.0) go to 20
72  gamit = 0.0
73  IF (ainta.GT.0.0 .OR. aeps.NE.0.0) gamit = gamr(a+1.0)
74  RETURN
75 C
76  20 IF (x.GT.1.0) go to 40
77  IF (a.GE.(-0.5) .OR. aeps.NE.0.0) CALL algams(a+1.0, algap1,
78  1 sgngam)
79  gamit = r9gmit(a, x, algap1, sgngam, alx)
80  RETURN
81 C
82  40 IF (a.LT.x) go to 50
83  t = r9lgit(a, x, alngam(a+1.0))
84  IF (t.LT.bot) CALL xerclr
85  gamit = exp(t)
86  RETURN
87 C
88  50 alng = r9lgic(a, x, alx)
89 C
90 C EVALUATE GAMIT IN TERMS OF LOG(GAMIC(A,X))
91 C
92  h = 1.0
93  IF (aeps.EQ.0.0 .AND. ainta.LE.0.0) go to 60
94  CALL algams(a+1.0, algap1, sgngam)
95  t = log(abs(a)) + alng - algap1
96  IF (t.GT.alneps) go to 70
97  IF (t.GT.(-alneps)) h = 1.0 - sga*sgngam*exp(t)
98  IF (abs(h).GT.sqeps) go to 60
99  CALL xerclr
100  CALL xermsg('SLATEC', 'GAMIT', 'RESULT LT HALF PRECISION', 1, 1)
101 C
102  60 t = -a*alx + log(abs(h))
103  IF (t.LT.bot) CALL xerclr
104  gamit = sign(exp(t), h)
105  RETURN
106 C
107  70 t = t - a*alx
108  IF (t.LT.bot) CALL xerclr
109  gamit = -sga*sgngam*exp(t)
110  RETURN
111 C
112  END