Blame view
fvn_fnlib/r9gmic.f
1.7 KB
38581db0c
git-svn-id: https...
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 |
function r9gmic (a, x, alx)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
c
c compute the complementary incomplete gamma function for a near
c a negative integer and for small x.
c
external alngam, r1mach
data euler / .5772156649 015329 e0 /
data eps, bot / 2*0.0 /
c
if (eps.eq.0.0) eps = 0.5*r1mach(3)
if (bot.eq.0.0) bot = alog(r1mach(1))
c
if (a.gt.0.0) call seteru (
1 41hr9gmic a must be near a negative integer, 41, 2, 2)
if (x.le.0.0) call seteru (25hr9gmic x must be gt zero, 25, 3, 2)
c
ma = a - 0.5
fm = -ma
m = -ma
c
te = 1.0
t = 1.0
s = t
do 20 k=1,200
fkp1 = k + 1
te = -x*te/(fm+fkp1)
t = te/fkp1
s = s + t
if (abs(t).lt.eps*s) go to 30
20 continue
call seteru (57hr9gmic no convergence in 200 terms of continued f
1raction, 57, 4, 2)
c
30 r9gmic = -alx - euler + x*s/(fm+1.0)
if (m.eq.0) return
c
if (m.eq.1) r9gmic = -r9gmic - 1.0 + 1.0/x
if (m.eq.1) return
c
te = fm
t = 1.0
s = t
mm1 = m - 1
do 40 k=1,mm1
fk = k
te = -x*te/fk
t = te/(fm-fk)
s = s + t
if (abs(t).lt.eps*abs(s)) go to 50
40 continue
c
50 do 60 k=1,m
r9gmic = r9gmic + 1.0/float(k)
60 continue
c
sgng = 1.0
if (mod(m,2).eq.1) sgng = -1.0
alng = alog(r9gmic) - alngam(fm+1.0)
c
r9gmic = 0.0
if (alng.gt.bot) r9gmic = sgng*exp(alng)
if (s.ne.0.0) r9gmic = r9gmic + sign (exp(-fm*alx+alog(abs(s)/fm))
1 , s)
c
if (r9gmic.eq.0.0 .and. s.eq.0.0) call seteru (
1 25hr9gmic result underflows, 25, 1, 0)
return
c
end
|