lxsd / fvn

Blame view

fvn_fnlib/zpsi.f 2.83 KB
edit raw normal view history
0c3098aed   cwaterkeyn   ChW 02/2010 for t... 
1
 
        complex(kind(1.d0)) function zpsi (zin)
38581db0c   daniau   git-svn-id: https... 
2
3
 
        implicit none
  c may 1978 edition.  w. fullerton, c3, los alamos scientific lab.
0c3098aed   cwaterkeyn   ChW 02/2010 for t... 
4
 
        complex(kind(1.d0)) zin, z, z2inv, corr,  zcot
38581db0c   daniau   git-svn-id: https... 
5
 
        dimension bern(13)
0c3098aed   cwaterkeyn   ChW 02/2010 for t... 
6
7
 
        real(kind(1.d0)) bern,d1mach,pi,bound, dxrel, rmin, rbig
        real(kind(1.d0)) x,y,cabsz
38581db0c   daniau   git-svn-id: https... 
8
9
10
 
        integer nterm,ndx,n,i
        external zcot, d1mach
  c
00055ac08   kwagner   Define some const... 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 
        data bern( 1) /   .8333333333 3333333 d-1 /
        data bern( 2) /  -.8333333333 3333333 d-2 /
        data bern( 3) /   .3968253968 2539683 d-2 /
        data bern( 4) /  -.4166666666 6666667 d-2 /
        data bern( 5) /   .7575757575 7575758 d-2 /
        data bern( 6) /  -.2109279609 2796093 d-1 /
        data bern( 7) /   .8333333333 3333333 d-1 /
        data bern( 8) /  -.4432598039 2156863 d0 /
        data bern( 9) /   .3053954330 2701197 d1 /
        data bern(10) /  -.2645621212 1212121 d2 /
        data bern(11) /   .2814601449 2753623 d3 /
        data bern(12) /  -.3454885393 7728938 d4 /
        data bern(13) /   .5482758333 3333333 d5 /
  c
        data pi / 3.141592653 589793 d0 /
        data nterm, bound, dxrel, rmin, rbig / 0d0, 4*0.0d0 /
38581db0c   daniau   git-svn-id: https... 
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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
 
  c
        if (nterm.ne.0) go to 10
        nterm = -0.30*log(d1mach(3))
  c maybe bound = n*(0.1*eps)**(-1/(2*n-1)) / (pi*exp(1))
        bound = 0.1171*dble(nterm) *
       1  (0.1*d1mach(3))**(-1.0/(2.0*dble(nterm)-1.0))
        dxrel = sqrt(d1mach(4))
        rmin = exp (max(log(d1mach(1)), -log(d1mach(2))) + 0.011 )
        rbig = 1.0/d1mach(3)
  c
   10   z = zin
        x = real(z)
        y = aimag(z)
        if (y.lt.0.0) z = conjg(z)
  c
        corr = (0.0, 0.0)
        cabsz = abs(z)
        if (x.ge.0.0 .and. cabsz.gt.bound) go to 50
        if (x.lt.0.0 .and. abs(y).gt.bound) go to 50
  c
        if (cabsz.lt.bound) go to 20
  c
  c use the reflection formula for real(z) negative, cabs(z) large, and
  c abs(aimag(y)) small.
  c
        corr = -pi*zcot(pi*z)
        z = 1.0 - z
        go to 50
  c
  c use the recursion relation for cabs(z) small.
  c
   20   if (cabsz.lt.rmin) call seteru (
       1  56hzpsi    zpsi called with z so near 0 that zpsi overflows,
       1  56, 2, 2)
  c
        if (x.ge.(-0.5) .or. abs(y).gt.dxrel) go to 30
        if ( abs((z-aint(x-0.5))/x).lt.dxrel) call seteru (
       1  68hzpsi    answer lt half precision because z too near negative
       2integer, 68, 1, 1)
        if (y.eq.0.0 .and. x.eq.aint(x)) call seteru (
       1  31hzpsi    z is a negative integer, 31, 3, 2)
  c
   30   n = sqrt(bound**2-y**2) - x + 1.0
        do 40 i=1,n
          corr = corr - 1.0/z
          z = z + 1.0
   40   continue
  c
  c now evaluate the asymptotic series for suitably large z.
  c
   50   if (cabsz.gt.rbig) zpsi = log(z) + corr
        if (cabsz.gt.rbig) go to 70
  c
        zpsi = (0.0, 0.0)
        z2inv = 1.0/z**2
        do 60 i=1,nterm
          ndx = nterm + 1 - i
          zpsi = bern(ndx) + z2inv*zpsi
   60   continue
        zpsi = log(z) - 0.5/z - zpsi*z2inv + corr
  c
   70   if (y.lt.0.0) zpsi = conjg(zpsi)
  c
        return
        end