Blame view

fvn_fnlib/zexprl.f 1.33 KB
0c3098aed   cwaterkeyn   ChW 02/2010 for t...
1
        complex(kind(1.d0)) function zexprl (z)
38581db0c   daniau   git-svn-id: https...
2
3
4
5
6
7
8
9
10
        implicit none
  c august 1980 edition.  w. fullerton, c3, los alamos scientific lab.
  c
  c evaluate  (cexp(z)-1)/z .  for small cabs(z), we use the taylor
  c series.  we could instead use the expression
  c        cexprl(z) = (exp(x)*exp(i*y)-1)/z
  c                  = (x*exprel(x) * (1 - 2*sin(y/2)**2) - 2*sin(y/2)**2
  c                                    + i*sin(y)*(1+x*exprel(x))) / z
  c
0c3098aed   cwaterkeyn   ChW 02/2010 for t...
11
12
        complex(kind(1.d0)) z
        real(kind(1.d0)) sqeps,r,xn,xln,rbnd,d1mach,alneps
38581db0c   daniau   git-svn-id: https...
13
14
        external d1mach
        integer nterms,irold,irold2,i
00055ac08   kwagner   Define some const...
15
        data nterms, rbnd, sqeps / 0d0, 2*0.0d0 /
38581db0c   daniau   git-svn-id: https...
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
  c
        if (nterms.ne.0) go to 10
        alneps = log(d1mach(3))
        xn = 3.72 - 0.3*alneps
        xln = log((xn+1.0)/1.36)
        nterms = xn - (xn*xln+alneps)/(xln+1.36) + 1.5
        rbnd = d1mach(3)
        sqeps = 0.0
  c
   10   r = abs(z)
        if (r.le.0.5) go to 20
  c
        call entsrc (irold, 1)
        zexprl = (exp(z) - 1.0) / z
        call erroff
        call entsrc (irold2, irold)
        if (abs(real(z)).lt.sqeps .and.  abs(zexprl).lt.sqeps) call
       1  seteru (41hzexprl  answer lt half prec, z near ni2pi,
       2  41, 1, 1)
        return
  c
   20   zexprl = (1.0, 0.0)
        if (r.lt.rbnd) return
  c
        zexprl = (0.0, 0.0)
        do 30 i=1,nterms
          zexprl = 1.0 + zexprl*z/dble(nterms+2-i)
   30   continue
  c
        return
        end