spbdphy.f 2.01 KB
 Sebastian Heimann committed Mar 31, 2015 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 `````` Sebastian Heimann committed Apr 01, 2015 11 ``````c dphy(n,xa,xb)=(phy(n,xa)-phy(n,xb))*2*(2*n-1)/(xa^2-xb^2) c `````` Sebastian Heimann committed Mar 31, 2015 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 `````` Sebastian Heimann committed Apr 01, 2015 21 22 `````` integer n double complex xa,xb `````` Sebastian Heimann committed Mar 31, 2015 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 `````` Sebastian Heimann committed Apr 01, 2015 36 37 38 `````` ca=(1.d0,0.d0) gn=ca fn=ca `````` Sebastian Heimann committed Mar 31, 2015 39 40 `````` cy=fn 100 i=i+1 `````` Sebastian Heimann committed Nov 24, 2016 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 committed Mar 31, 2015 44 45 `````` cy=cy+fn if(cdabs(fn).gt.eps*cdabs(cy))goto 100 `````` Sebastian Heimann committed Apr 01, 2015 46 `````` spbdphy=cy `````` Sebastian Heimann committed Mar 31, 2015 47 `````` return `````` Sebastian Heimann committed Jul 19, 2013 48 `` end``