lxsd / fvn

Blame view

fvn_fnlib/d9gmic.f 1.88 KB
edit raw normal view history
38581db0c   daniau   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
71
72
 
        double precision function d9gmic (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
        double precision a, x, alx, alng, bot, eps, euler, fk, fkp1, fm,
       1  s, sgng, t, te, d1mach, dlngam, dexp, dlog
        external d1mach, dlngam
        data euler / 0.5772156649 0153286060 6512090082 40 d0 /
        data eps, bot / 2*0.d0 /
  c
        if (eps.ne.0.d0) go to 10
        eps = 0.5d0*d1mach(3)
        bot = dlog (d1mach(1))
  c
   10   if (a.gt.0.d0) call seteru (
       1  41hd9gmic  a must be near a negative integer, 41, 2, 2)
        if (x.le.0.d0) call seteru (25hd9gmic  x must be gt zero, 25, 3,2)
  c
        m = -(a - 0.5d0)
        fm = m
  c
        te = 1.0d0
        t = 1.0d0
        s = t
        do 20 k=1,200
          fkp1 = k + 1
          te = -x*te/(fm+fkp1)
          t = te/fkp1
          s = s + t
          if (dabs(t).lt.eps*s) go to 30
   20   continue
        call seteru (57hd9gmic  no convergence in 200 terms of continued f
       1raction, 57, 4, 2)
  c
   30   d9gmic = -alx - euler + x*s/(fm+1.0d0)
        if (m.eq.0) return
  c
        if (m.eq.1) d9gmic = -d9gmic - 1.d0 + 1.d0/x
        if (m.eq.1) return
  c
        te = fm
        t = 1.d0
        s = t
        mm1 = m - 1
        do 40 k=1,mm1
          fk = k
          te = -x*te/fk
          t = te/(fm-fk)
          s = s + t
          if (dabs(t).lt.eps*dabs(s)) go to 50
   40   continue
  c
   50   do 60 k=1,m
          d9gmic = d9gmic + 1.0d0/dble(float(k))
   60   continue
  c
        sgng = 1.0d0
        if (mod(m,2).eq.1) sgng = -1.0d0
        alng = dlog(d9gmic) - dlngam(fm+1.d0)
  c
        d9gmic = 0.d0
        if (alng.gt.bot) d9gmic = sgng * dexp(alng)
        if (s.ne.0.d0) d9gmic = d9gmic +
       1  dsign (dexp(-fm*alx+dlog(dabs(s)/fm)), s)
  c
        if (d9gmic.eq.0.d0 .and. s.eq.0.d0) call seteru (
       1  25hd9gmic  result underflows, 25, 1, 0)
        return
  c
        end