1!  *************************************************************************
2!
3!        This is a slightly modified version of Log(Gamma) function
4!        program from Num. Rec.
5!       Serves as backup if intrinic Gamma function
6!       is not available
7!  ************************************************************************
8
9      DOUBLE PRECISION FUNCTION PAW_LN_GAMMA (XX)
10      implicit none
11      DOUBLE PRECISION XX
12      INTEGER J
13      DOUBLE PRECISION SER
14      DOUBLE PRECISION STP
15      DOUBLE PRECISION TMP
16      DOUBLE PRECISION X
17      DOUBLE PRECISION Y, COF(6)
18      SAVE STP, COF
19      DATA COF, STP/ 76.18009172947146D0, -86.50532032941677D0,
20     1   24.01409824083091D0, -1.231739572450155D0,
21     2   0.1208650973866179D-2, -.5395239384953D-5, 2.5066282746310005D0
22     3   /
23      X = XX
24      Y = X
25      TMP = X + 5.5D0
26      TMP = (X + 0.5D0)*DLOG(TMP) - TMP
27      SER = 1.000000000190015D0
28      DO J = 1, 6
29         Y = Y + 1.0D0
30         SER = SER + COF(J)/Y
31      END DO
32
33      PAW_LN_GAMMA = TMP + DLOG(STP*SER/X)
34
35      RETURN
36      END
37
38!  *************************************************
39!
40!     Name    :
41!
42!
43!     Purpose :
44!
45!
46!     Created :
47!
48!  *************************************************
49      DOUBLE PRECISION FUNCTION PAW_GAMMA(X)
50      implicit none
51      DOUBLE PRECISION X
52
53      DOUBLE PRECISION XX
54      DOUBLE PRECISION PAW_LN_GAMMA
55
56      XX = X
57      PAW_GAMMA = DEXP(PAW_LN_GAMMA(XX))
58
59      return
60      END
61
62
63
64!  *************************************************
65!
66!     Name    :
67!
68!
69!     Purpose :
70!
71!
72!     Created :
73!
74!  *************************************************
75      DOUBLE PRECISION FUNCTION PAW_GAMMP (A, X)
76      implicit none
77      DOUBLE PRECISION A, X
78
79      DOUBLE PRECISION GAMMCF, GAMSER, GLN
80
81      IF (X .LT. A+1.0D0) THEN
82
83         CALL PAW_GSER(GAMSER, A, X, GLN)
84         PAW_GAMMP = GAMSER
85
86      ELSE
87
88         CALL PAW_GCF(GAMMCF, A, X, GLN)
89         PAW_GAMMP = 1.0D0 - GAMMCF
90
91      ENDIF
92
93      return
94      END
95
96!  *************************************************
97!
98!     Name    :
99!
100!
101!     Purpose :
102!
103!
104!     Created :
105!
106!  *************************************************
107      SUBROUTINE PAW_GCF(GAMMCF, A, X, GLN)
108      implicit none
109      INTEGER ITMAX
110      DOUBLE PRECISION A, GAMMCF, GLN, X, EPS, FPMIN
111      PARAMETER (ITMAX = 100, EPS = 3.D-16, FPMIN = 1.D-30)
112      DOUBLE PRECISION AN, B, C, D, DEL, H
113      INTEGER I
114
115      DOUBLE PRECISION PAW_LN_GAMMA
116      EXTERNAL         PAW_LN_GAMMA
117
118      GLN = PAW_LN_GAMMA(A)
119      B = X + 1.0D0 - A
120      C = 1.0D0/1.D-30
121      D = 1.0D0/B
122      H = D
123      DO I = 1, 100
124         AN = -I*(I - A)
125         B = B + 2.0D0
126         D = AN*D + B
127         IF (DABS(D) .LT. 1.D-30) D = 1.D-30
128         C = B + AN/C
129         IF (DABS(C) .LT. 1.D-30) C = 1.D-30
130         D = 1.0D0/D
131         DEL = D*C
132         H = H*DEL
133         IF (DABS(DEL - 1.0D0) .LT. 3.D-16) GO TO 1
134      END DO
135      PAUSE 'a too large, ITMAX too small in gcf'
136    1 CONTINUE
137      GAMMCF = DEXP((-X) + A*DLOG(X) - GLN)*H
138
139      return
140      END
141
142!  *************************************************
143!
144!     Name    :
145!
146!
147!     Purpose :
148!
149!
150!     Created :
151!
152!  *************************************************
153      SUBROUTINE PAW_GSER(GAMSER, A, X, GLN)
154      implicit none
155      DOUBLE PRECISION A, X
156      DOUBLE PRECISION GAMSER, GLN
157
158!    *** local variables ***
159      INTEGER ITMAX
160      PARAMETER (ITMAX = 100)
161      DOUBLE PRECISION EPS
162      PARAMETER (EPS = 3.0D-16)
163      INTEGER N
164      DOUBLE PRECISION AP, DEL, SUM
165
166      DOUBLE PRECISION PAW_LN_GAMMA
167      EXTERNAL         PAW_LN_GAMMA
168
169
170      GLN = PAW_LN_GAMMA(A)
171
172      IF (X .LE. 0.0D0) THEN
173      IF(X.lt.0.0d0) CALL errquit("x < 0 in PAW_GSER",0,1)
174         GAMSER = 0.0D0
175         RETURN
176      ENDIF
177
178      AP = A
179      SUM = 1.0D0/A
180      DEL = SUM
181      DO N = 1, 100
182         AP = AP + 1.0D0
183         DEL = DEL*X/AP
184         SUM = SUM + DEL
185         IF (DABS(DEL) .LT. DABS(SUM)*3.0D-16) GO TO 1
186
187      END DO
188
189      CALL errquit
190     >     ("a too large,ITMAX too small in PAW_GSER",0,1)
191
192    1 CONTINUE
193      GAMSER = SUM*DEXP((-X) + A*DLOG(X) - GLN)
194
195      return
196      END
197
198
199
200c $Id$
201