spbdphy.f 2.01 KB
Newer Older
Sebastian Heimann's avatar
Sebastian Heimann committed
1
2
3
4
5
6
7
8
9
10
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Function: spbdpsy                                                         c
c     calculates eq.(103) of Takeuchi & Saito related to spherical Bessel       c
c     functions for complex argument                                            c
c                                                                               c
c     Input:                                                                    c
c     n = order                                                                 c
c     x = complex argument                                                      c
c                                                                               c
c     Return:                                                                   c
11
c     dphy(n,xa,xb)=(phy(n,xa)-phy(n,xb))*2*(2*n-1)/(xa^2-xb^2)                 c
Sebastian Heimann's avatar
Sebastian Heimann committed
12
13
14
15
16
17
18
19
20
c     for |x/2|^2 << n                                                          c
c                                                                               c
c     First implemention: 17 June 2007                                          c
c     Rongjiang Wang                                                            c
c     GeoForschungsZentrum Potsdam, Germany                                     c
c     Email: wang@gfz-potsdam.de                                                c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double complex function spbdphy(n,xa,xb)
      implicit none
21
22
      integer n
      double complex xa,xb
Sebastian Heimann's avatar
Sebastian Heimann committed
23
24
25
26
27
28
29
30
31
32
33
34
35
c
c     local memory
c
      integer i
      double complex cxa2,cxb2,fn,gn,ca,cy
c
      double precision eps
      data eps/1.0d-08/
c
      cxa2=(0.5d0,0.d0)*xa*xa
      cxb2=(0.5d0,0.d0)*xb*xb
c
      i=1
36
37
38
      ca=(1.d0,0.d0)
      gn=ca
      fn=ca
Sebastian Heimann's avatar
Sebastian Heimann committed
39
40
      cy=fn
100   i=i+1
41
42
43
      ca=(1.d0,0.d0)/dcmplx(dble((2*(n-i)+1)*i),0.d0)
      gn=gn*cxb2*ca
      fn=gn+fn*cxa2*ca
Sebastian Heimann's avatar
Sebastian Heimann committed
44
45
      cy=cy+fn
      if(cdabs(fn).gt.eps*cdabs(cy))goto 100
46
      spbdphy=cy
Sebastian Heimann's avatar
Sebastian Heimann committed
47
      return
Sebastian Heimann's avatar
Sebastian Heimann committed
48
      end