lxsd / fvn

Blame view

fvn_fnlib/d9chu.f 2.11 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 d9chu (a, b, z)
  c august 1977 edition.  w. fullerton, c3, los alamos scientific lab.
  c
  c evaluate for large z  z**a * u(a,b,z)  where u is the logarithmic
  c confluent hypergeometric function.  a rational approximation due to y.
  c l. luke is used.  when u is not in the asymptotic region, i.e., when a
  c or b is large compared with z, considerable significance loss occurs.
  c a warning is provided when the computed result is less than half
  c precision.
  c
        double precision a, b, z, aa(4), bb(4), ab, anbn, bp, ct1, ct2,
       1  ct3, c2, d1z, eps, g1, g2, g3, sab, sqeps, x2i1,  d1mach,
       2  dsqrt
        external d1mach
        data eps, sqeps / 2*0.0d0 /
  c
        if (eps.ne.0.0d0) go to 10
        eps = 4.0d0*d1mach(4)
        sqeps = dsqrt (d1mach(4))
  c
   10   bp = 1.0d0 + a - b
        ab = a*bp
        ct2 = 2.0d0 * (z - ab)
        sab = a + bp
  c
        bb(1) = 1.0d0
        aa(1) = 1.0d0
  c
        ct3 = sab + 1.0d0 + ab
        bb(2) = 1.0d0 + 2.0d0*z/ct3
        aa(2) = 1.0d0 + ct2/ct3
  c
        anbn = ct3 + sab + 3.0d0
        ct1 = 1.0d0 + 2.0d0*z/anbn
        bb(3) = 1.0d0 + 6.0d0*ct1*z/ct3
        aa(3) = 1.0d0 + 6.0d0*ab/anbn + 3.0d0*ct1*ct2/ct3
  c
        do 30 i=4,300
          x2i1 = 2*i - 3
          ct1 = x2i1/(x2i1-2.0d0)
          anbn = anbn + x2i1 + sab
          ct2 = (x2i1 - 1.0d0)/anbn
          c2 = x2i1*ct2 - 1.0d0
          d1z = x2i1*2.0d0*z/anbn
  c
          ct3 = sab*ct2
          g1 = d1z + ct1*(c2+ct3)
          g2 = d1z - c2
          g3 = ct1*(1.0d0 - ct3 - 2.0d0*ct2)
  c
          bb(4) = g1*bb(3) + g2*bb(2) + g3*bb(1)
          aa(4) = g1*aa(3) + g2*aa(2) + g3*aa(1)
  c
          d9chu = aa(4)/bb(4)
          if (dabs(d9chu-aa(1)/bb(1)).lt.eps*dabs(d9chu)) go to 40
  c
  c if overflows or underflows prove to be a problem, the statements
  c below could be altered to incorporate a dynamically adjusted scale
  c factor.
  c
          do 20 j=1,3
            aa(j) = aa(j+1)
            bb(j) = bb(j+1)
   20     continue
   30   continue
        call seteru (35hd9chu   no convergence in 300 terms, 35, 1, 2)
  c
   40     if (d9chu.lt.sqeps .or. d9chu.gt.1.0/sqeps) call seteru (
       1  32hd9chu   answer lt half precision, 32, 2, 1)
  c
        return
        end