1*     ******************************************************
2*     *                                                    *
3*     *             nwpw_SpecialKummer                     *
4*     *                                                    *
5*     ******************************************************
6*
7*     Calculates a special case of the Kummer confluent hypergeometric
8*     function, M(n+1/2,l+3/2,z) for z .LE. 0
9*
10*     This function was created by  Marat Valiev, and  modified by Eric Bylaska.
11*     See Abramowitz and Stegun for the formulas used in this function.
12*
13      real*8 function nwpw_SpecialKummer(n,l,z)
14      implicit none
15      integer n,l
16      real*8  z
17
18*     *** local variables ***
19      real*8 eps
20      parameter (eps=1.0d-16)
21
22      integer i
23      real*8 a,b,m1,m3,s
24
25*     **** external functions ****
26      real*8   nwpw_gamma,nwpw_gammp
27      external nwpw_gamma,nwpw_gammp
28
29      nwpw_SpecialKummer = 0.0d0
30
31*    *** cannot handle positive z ***
32      if (z.gt.0.0d0) then
33         call errquit('nwpw_SpecialKummer:invalid parameter, z>0',0,1)
34      end if
35
36
37*    *** solution for z==0 ***
38      if (z.eq.0.0d0) then
39         nwpw_SpecialKummer = 1.0d0
40         return
41      end if
42
43*     ***** M(a,a+1,z) = a * (-z)**(-a) * igamma(a,-z) = a * (-z)**(-a) * P(a,-z) *Gamma(a)  where z is real and a = (n+0.5)  ****
44      if (n.eq.l) then
45         nwpw_SpecialKummer = nwpw_gammp(n+0.5d0,(-z))
46     >                       *(n+0.5d0)
47     >                       *((-z)**((-n)- 0.5d0))
48     >                       *nwpw_gamma(n+0.5d0)
49         return
50
51*     ***** M(a,a,z) = exp(z)  where a = (n+0.5)  ****
52      else if (n.eq.(l+1)) then
53         nwpw_SpecialKummer = dexp(z)
54         return
55      end if
56
57!     *** do inifinite series for small z
58      if (dabs(z).le.1.0d0) then
59
60         nwpw_SpecialKummer = 1.0d0
61         s = 1.0d0
62         a = n + 0.5d0
63         b = l + 1.5d0
64         do i=1,10000
65            s = s*(a+i-1)*z/((b+i-1)*i)
66            nwpw_SpecialKummer = nwpw_SpecialKummer + s
67            if (dabs(s).lt.eps) return
68         end do
69         call errquit("nwpw_SpecialKummer:cannot converge",0,1)
70         return
71      end if
72
73      if (n.lt.l) then
74
75      !*** starting point n=l or b=a+1***
76         a = n + 0.5d0
77         b = n + 1.5d0
78
79      !*** m1 = M(a,b-1) ***
80      !*** m2 = M(a,b,z) ***
81         m1 = dexp(z)
82         nwpw_SpecialKummer = nwpw_gammp(a,(-z))*a/(-z)**a*nwpw_gamma(a)
83
84      !**********************************************
85      ! using recursion formula
86      ! z(a-b)M(a,b+1,z)=b(b-1)M(a,b-1,z)+b(1-b-z)M(a,b,z)
87      ! obtain M(1/2,3/2+l  ,z) --> m2
88      !        M(1/2,3/2+l-1,z) --> m2
89      !**********************************************
90         do i=1,l-n
91            m3=(b*(b-1.0d0)*m1+b*(1.0d0-b-z)*nwpw_SpecialKummer)
92     >         /(z*(a-b))
93            b = b + 1
94            m1 = nwpw_SpecialKummer
95            nwpw_SpecialKummer = m3
96         end do
97
98      else if (n.gt.(l+1)) then
99
100      !*** starting point n=l+1 or b=a ***
101         a = l + 1.5d0
102         b = l + 1.5d0
103
104      !*** m1 = M(a-1,b) ***
105      !*** m2 = M(a,a,z) ***
106         m1 = nwpw_gammp(a-1.0d0,(-z))*(a-1.0d0)/(-z)**(a-1.0d0)*
107     >      nwpw_gamma(a-1.0d0)
108         nwpw_SpecialKummer = dexp(z)
109
110      !**********************************************
111      ! using recursion formula
112      ! aM(a+1,b,z)=(b-a)M(a-1,b,z)+(2a-b+z)M(a,b,z)
113      ! obtain M(n+1/2-1,3/2,z)   --> m1
114      !        M(n+1/2  ,3/2,z)   --> m2
115      !**********************************************
116         do i=1,n-l-1
117            m3 = ((b-a)*m1+(2*a-b+z)*nwpw_SpecialKummer)/a
118            m1 = nwpw_SpecialKummer
119            nwpw_SpecialKummer = m3
120            a = a + 1
121         end do
122      end if
123
124      return
125      end
126
127*     ***************************************
128*     *                                     *
129*     *          nwpw_ln_gamma              *
130*     *                                     *
131*     ***************************************
132
133*        This function computes the Log(Gamma) function
134*       Serves as backup if intrinic Gamma function is not available
135
136      real*8 function nwpw_ln_gamma(XX)
137      implicit none
138      real*8 XX
139      integer J
140      real*8 SER
141      real*8 STP
142      real*8 TMP
143      real*8 X
144      DOUBLE PRECISION Y, COF(6)
145      SAVE STP, COF
146      DATA COF, STP/ 76.18009172947146d0, -86.50532032941677d0,
147     1   24.01409824083091d0, -1.231739572450155d0,
148     2   0.1208650973866179d-2, -.5395239384953d-5, 2.5066282746310005d0
149     3   /
150      X = XX
151      Y = X
152      TMP = X + 5.5d0
153      TMP = (X + 0.5d0)*dlog(TMP) - TMP
154      SER = 1.000000000190015d0
155      DO J = 1, 6
156         Y = Y + 1.0d0
157         SER = SER + COF(J)/Y
158      END DO
159      nwpw_ln_gamma = TMP + dlog(STP*SER/X)
160      return
161      end
162
163*     **********************************************
164*     *                                            *
165*     *            nwpw_gamma                      *
166*     *                                            *
167*     **********************************************
168*
169*     This function computes the Gamma function
170*
171      real*8 function nwpw_gamma(X)
172      implicit none
173      real*8 X
174
175      real*8 XX
176      real*8 nwpw_ln_gamma
177      external nwpw_ln_gamma
178
179      XX = X
180      nwpw_gamma = dexp(nwpw_ln_gamma(XX))
181      return
182      end
183
184*     **********************************************
185*     *                                            *
186*     *           nwpw_gammp                       *
187*     *                                            *
188*     **********************************************
189*
190*     This function computes the incomplete gamma function P(a,x) = igamma(a,x)/Gamma(a)
191*
192*                             /x
193*       P(a,x) = 1/Gamma(a) * | exp(-t) * t**(a-1) dt
194*                             /0
195*
196      real*8 function nwpw_gammp(A, X)
197      implicit none
198      real*8 a,x
199      real*8 gammcf,gamser,gln
200
201      if (x.lt. (a+1.0d0)) then
202         call nwpw_gser(gamser, a, x, gln)
203         nwpw_gammp = gamser
204      else
205         call nwpw_gcf(gammcf,a,x,gln)
206         nwpw_gammp = 1.0d0 - gammcf
207      end if
208      return
209      end
210
211*     **********************************************
212*     *                                            *
213*     *               nwpw_gcf                     *
214*     *                                            *
215*     **********************************************
216      subroutine nwpw_gcf(GAMMCF, A, X, GLN)
217      implicit none
218      integer ITMAX
219      real*8 A,GAMMCF,GLN,X,EPS,FPMIN
220      parameter (ITMAX = 100, EPS = 3.0d-16, FPMIN = 1.0d-30)
221      real*8 AN,B,C,D,DEL,H
222      integer I
223
224      real*8   nwpw_ln_gamma
225      external nwpw_ln_gamma
226
227      GLN = nwpw_ln_gamma(A)
228      B = X + 1.0d0 - A
229      C = 1.0d0/1.0d-30
230      D = 1.0d0/B
231      H = D
232      DO I = 1, 100
233         AN = -I*(I - A)
234         B = B + 2.0d0
235         D = AN*D + B
236         IF (dabs(D) .LT. 1.0d-30) D = 1.0d-30
237         C = B + AN/C
238         IF (dabs(C) .LT. 1.0d-30) C = 1.0d-30
239         D = 1.0d0/D
240         DEL = D*C
241         H = H*DEL
242         IF (dabs(DEL-1.0d0) .LT. 3.0d-16) GO TO 1
243      END DO
244      call errquit('a too large, ITMAX too small in nwpw_gcf',0,0)
245
246    1 continue
247      GAMMCF = dexp((-X) + A*dlog(X) - GLN)*H
248      return
249      end
250
251*     **********************************************
252*     *                                            *
253*     *             nwpw_gser                      *
254*     *                                            *
255*     **********************************************
256      subroutine nwpw_gser(GAMSER, A, X, GLN)
257      implicit none
258      real*8 A,X
259      real*8 GAMSER,GLN
260
261!    *** local variables ***
262      integer ITMAX
263      parameter (ITMAX = 100)
264      real*8 EPS
265      parameter (EPS = 3.0d-16)
266      integer N
267      real*8 AP,DEL,SUM
268
269      real*8   nwpw_ln_gamma
270      external nwpw_ln_gamma
271
272      GLN = nwpw_ln_gamma(A)
273
274      if (X .le. 0.0d0) then
275         if (X.lt.0.0d0) call errquit("x < 0 in nwpw_gser",0,1)
276         GAMSER = 0.0d0
277         return
278      end if
279
280      AP = A
281      SUM = 1.0d0/A
282      DEL = SUM
283      do N = 1, 100
284         AP = AP + 1.0d0
285         DEL = DEL*X/AP
286         SUM = SUM + DEL
287         if (dabs(DEL) .lt. dabs(SUM)*3.0d-16) GO TO 1
288
289      end do
290
291      call errquit
292     >     ("a too large,ITMAX too small in nwpw_gser",0,1)
293
294    1 continue
295      GAMSER = SUM*dexp((-X) + A*dlog(X) - GLN)
296
297      return
298      end
299
300c $Id$
301