lxsd / fvn

Blame view

fvn_fnlib/r9gmit.f 1.46 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
 
        function r9gmit (a, x, algap1, sgngam, alx)
  c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
  c
  c compute tricomi-s incomplete gamma function for small x.
  c
        external alngam, r1mach
        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 (x.le.0.0) call seteru (24hr9gmit  x should be gt 0, 24, 1, 2)
  c
        ma = a + 0.5
        if (a.lt.0.0) ma = a - 0.5
        aeps = a - float(ma)
  c
        ae = a
        if (a.lt.(-0.5)) ae = aeps
  c
        t = 1.0
        te = ae
        s = t
        do 20 k=1,200
          fk = k
          te = -x*te/fk
          t = te/(ae+fk)
          s = s + t
          if (abs(t).lt.eps*abs(s)) go to 30
   20   continue
        call seteru (54hr9gmit  no convergence in 200 terms of taylor-s se
       1ries, 54, 2, 2)
  c
   30   if (a.ge.(-0.5)) algs = -algap1 + alog(s)
        if (a.ge.(-0.5)) go to 60
  c
        algs = -alngam(1.0+aeps) + alog(s)
        s = 1.0
        m = -ma - 1
        if (m.eq.0) go to 50
        t = 1.0
        do 40 k=1,m
          t = x*t/(aeps-float(m+1-k))
          s = s + t
          if (abs(t).lt.eps*abs(s)) go to 50
   40   continue
  c
   50   r9gmit = 0.0
        algs = -float(ma)*alog(x) + algs
        if (s.eq.0.0 .or. aeps.eq.0.0) go to 60
  c
        sgng2 = sgngam*sign(1.0,s)
        alg2 = -x - algap1 + alog(abs(s))
  c
        if (alg2.gt.bot) r9gmit = sgng2*exp(alg2)
        if (algs.gt.bot) r9gmit = r9gmit + exp(algs)
        return
  c
   60   r9gmit = exp(algs)
        return
  c
        end