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