C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Bessel and Hankel function routines, from Anders Bostrom C 891108, a few minor changes by Ilkka Karasalo (IK) C also, a few minor changes by Sven Ivansson (SI) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C THE FILE CONTAINS THE MAIN SUBROUTINES C C CDHAN C CDBES C DHAN C DBES C C WHICH COMPUTE HANKEL AND BESSEL FUNCTIONS C C CSUMNEUSI (adapted from csumneu) HAS BEEN ADDED BY SVEN IVANSSON C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C 891212: The functions CDBES and CDHAN now return zero C if ABS(Im(Z)) > ZIMMAX = 200. (This number was C found by numerical experiment.) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CDHAN C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE CDHAN(MMAX1,Z,B,H) C C COMPUTES THE BESSEL AND HANKEL FUNTIONS OF ORDER 0,1,...,MMAX1-1 C AND COMPLEX ARGUMENT Z AND PUTS THEM IN THE ARRAYS B(1),... AND H(1),... C MMAX1<52; AT LEAST 9 FIGURES ARE ACCURATE IN H AND 13 IN B C USES SUBROUTINES CDBES, CSUMBES, CASBES, CSUMNEU, AND CASHAN C implicit real*8 (a-h,o-z) complex*16 z,b(2),h(2),j0,j1,y0,y1,y2,h0,h1, * twobyz parameter ( zimmax = 200.d0) ! Max ABS(Im(Z)) allowed c ! See comment first in file logical*2 ovfl zab = abs(z) y=dimag(z) ovfl = abs(y).gt.zimmax if(ovfl.or.zab.eq.0.d0) then c Z = 0: Return H = 0 (H has pole at 0) b(1) = dcmplx(1.d0,0.d0) if(ovfl) b(1) = dcmplx(0.d0,0.d0) IF (ZAB.EQ.0.0D0) THEN C "overflow" but choose symmetry in " J0 = 0.5d0 * ( H0(1) + H0(2) ) " H(1) = (1.0D0,0.0D0) ELSE h(1) = dcmplx(0.d0,0.d0) END IF !!! SI do 1 i = 2, mmax1 b(i) = dcmplx(0.d0,0.d0) 1 h(i) = dcmplx(0.d0,0.d0) return end if x=dreal(z) twobyz = 2.0d0/z call cdbes(mmax1,z,b) if((x/12)**2+(y/8)**2.lt.1.d0) then j0=b(1) c If MMAX1 = 1, use a dummy value for J1. Note that c the Y0 value returned by CSUMNEU does not depend on c the input J1, and is correct if J0 is correct if(mmax1.ge.2) then j1=b(2) else j1 = dcmplx(0.d0,0.d0) end if call csumneu(z,j0,j1,y0,y1) h(1) = j0 + dcmplx(-dimag(y0),dreal(y0)) if(mmax1.ge.2) * h(2) = j1 + dcmplx(-dimag(y1),dreal(y1)) do 10 m=3,mmax1 y2 = twobyz*(m-2)*y1 - y0 y0 = y1 y1 = y2 10 h(m) = b(m) + dcmplx(-dimag(y1),dreal(y1)) else call cashan(z,h0,h1) h(1)=h0 if(mmax1.ge.2) h(2)=h1 do 30 m=3,mmax1 30 h(m)=twobyz*(m-2)*h(m-1)-h(m-2) end if 999 return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CDHANSVEN as cdhan but without "B" (only "H") C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE CDHANSVEN(MMAX1,Z,H) C C COMPUTES THE HANKEL FUNTIONS OF ORDER 0,1,...,MMAX1-1 C AND COMPLEX ARGUMENT Z AND PUTS THEM IN THE ARRAY H(1),... C MMAX1<52; AT LEAST 9 FIGURES ARE ACCURATE IN H C USES SUBROUTINES CDBES, CSUMBES, CASBES, CSUMNEU, AND CASHAN C implicit real*8 (a-h,o-z) complex*16 z,b(2),h(2),j0,j1,y0,y1,y2,h0,h1, * twobyz parameter ( zimmax = 200.d0) ! Max ABS(Im(Z)) allowed c ! See comment first in file logical*2 ovfl zab = abs(z) y=dimag(z) ovfl = y.lt.-zimmax if(ovfl.or.zab.eq.0.d0) then c Z = 0: Return H = 0 (H has pole at 0) IF (ZAB.EQ.0.0D0) THEN C "overflow" but choose symmetry in " J0 = 0.5d0 * ( H0(1) + H0(2) ) " H(1) = (1.0D0,0.0D0) ELSE h(1) = dcmplx(0.d0,0.d0) END IF !!! SI do 1 i = 2, mmax1 1 h(i) = dcmplx(0.d0,0.d0) return end if x=dreal(z) twobyz = 2.0d0/z if((x/12)**2+(y/8)**2.lt.1.d0) then call cdbes(mmax1,z,b) j0=b(1) c If MMAX1 = 1, use a dummy value for J1. Note that c the Y0 value returned by CSUMNEU does not depend on c the input J1, and is correct if J0 is correct if(mmax1.ge.2) then j1=b(2) else j1 = dcmplx(0.d0,0.d0) end if call csumneu(z,j0,j1,y0,y1) h(1) = j0 + dcmplx(-dimag(y0),dreal(y0)) if(mmax1.ge.2) * h(2) = j1 + dcmplx(-dimag(y1),dreal(y1)) do 10 m=3,mmax1 y2 = twobyz*(m-2)*y1 - y0 y0 = y1 y1 = y2 10 h(m) = b(m) + dcmplx(-dimag(y1),dreal(y1)) else call cashan(z,h0,h1) h(1)=h0 if(mmax1.ge.2) h(2)=h1 do 30 m=3,mmax1 30 h(m)=twobyz*(m-2)*h(m-1)-h(m-2) end if 999 return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CDHANSPEC as cdhansven but with "common /inspec/" active C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE CDHANSPEC(MMAX1,Z,H) C C COMPUTES THE "inspec scaled" HANKEL FUNTIONS OF ORDER 0,1,...,MMAX1-1 C AND COMPLEX ARGUMENT Z AND PUTS THEM IN THE ARRAY H(1),... C MMAX1<52; AT LEAST 9 FIGURES ARE ACCURATE IN H C USES SUBROUTINES CDBES, CSUMBES, CASBES, CSUMNEU, AND CASHAN C implicit real*8 (a-h,o-z) complex*16 z,b(2),h(2),j0,j1,y0,y1,y2,h0,h1, * twobyz INCLUDE '../cpm/INSPEC.INC' parameter ( zimmax = 200.d0) ! Max ABS(Im(Z)) allowed c ! See comment first in file logical*2 ovfl zab = abs(z) y=dimag(z) ovfl = y.lt.-zimmax if(ovfl.or.zab.eq.0.d0) then c Z = 0: Return H = 0 (H has pole at 0) IF (ZAB.EQ.0.0D0) THEN C "overflow" but choose symmetry in " J0 = 0.5d0 * ( H0(1) + H0(2) ) " IF (IRKSPEC.EQ.0) THEN H(1) = (1.0D0,0.0D0) ELSE H(1) = DCMPLX(DEXP(DLNSPEC),0.0D0) END IF ELSE h(1) = dcmplx(0.d0,0.d0) END IF !!! SI do 1 i = 2, mmax1 1 h(i) = dcmplx(0.d0,0.d0) return end if x=dreal(z) twobyz = 2.0d0/z if((x/12)**2+(y/8)**2.lt.1.d0) then call cdbes(mmax1,z,b) j0=b(1) c If MMAX1 = 1, use a dummy value for J1. Note that c the Y0 value returned by CSUMNEU does not depend on c the input J1, and is correct if J0 is correct if(mmax1.ge.2) then j1=b(2) else j1 = dcmplx(0.d0,0.d0) end if call csumneu(z,j0,j1,y0,y1) h(1) = j0 + dcmplx(-dimag(y0),dreal(y0)) if(mmax1.ge.2) * h(2) = j1 + dcmplx(-dimag(y1),dreal(y1)) do 10 m=3,mmax1 y2 = twobyz*(m-2)*y1 - y0 y0 = y1 y1 = y2 10 h(m) = b(m) + dcmplx(-dimag(y1),dreal(y1)) IF (IRKSPEC.NE.0) THEN SLASK = DEXP(DLNSPEC) DO M=1,MMAX1 H(M) = SLASK*H(M) END DO END IF else call cashanspec(z,h0,h1) h(1)=h0 if(mmax1.ge.2) h(2)=h1 do 30 m=3,mmax1 30 h(m)=twobyz*(m-2)*h(m-1)-h(m-2) end if 999 return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CASHAN C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine cashan(z,h0,h1) C C COMPUTES H0 AND H1 USING THE ASYMPTOTIC FORMULA C implicit real*8 (a-h,o-z) dimension ap0(13),ap1(13),aq0(13),aq1(13) complex*16 z,h0,h1,p0,p1,q0,q1,z2,a, onebyz, wor parameter ( twobypi = .6366197723675814D+00 ) ! 2/PI parameter ( pibyfour = .7853981633974483D+00 ) ! PI/4 parameter ( thrpibyf = .2356194490192345D+01 ) ! 3*PI/4 DATA AP0 /.1D1,-.703125D-1,.112152099609375D0, ,-.5725014209747314D0,.6074042001273483D1, ,-.1100171402692467D3,.3038090510922384D4, ,-.1188384262567832D6,.6252951493434794D7, ,-.4259392165047667D9,.3646840080706554D11, ,-.3833534661393942D13,.4854014686852896D15/ DATA AP1 /.1D1,.1171875D0,-.144195556640625D0, ,.6765925884246826D0,-.6883914268109948D1, ,.1215978918765359D3,-.3302272294480852D4, ,.1276412726461746D6,-.6656367718817685D7, ,.4502786003050390D9,-.3833857520742787D11, ,.4011838599133194D13,-.5060568503314721D15/ DATA AQ0 /-.125D0,.732421875D-1,-.2271080017089844D0, ,.1727727502584457D1,-.2438052969955606D2, ,.5513358961220205D3,-.1825775547429317D5, ,.8328593040162890D6,-.5006958953198890D8, ,.3836255180230431D10,-.3649010818849831D12, ,.4218971570284093D14,-.5827244631566901D16/ DATA AQ1 /.375D0,-.1025390625D0,.2775764465332031D0, ,-.1993531733751297D1,.2724882731126854D2, ,-.6038440767050701D3,.1971837591223663D5, ,-.8902978767070675D6,.5310411010968520D8, ,-.4043620325107752D10,.3827011346598603D12, ,-.4406481417852275D14,.6065091351222692D16/ IF (DIMAG(Z).LT.0.0D0) THEN IF (DREAL(Z).LT.0.0D0) THEN WRITE(*,*) '%%% besbod.f, better accuracy by symmetry; Sven I stop' STOP 5760 END IF END IF zab = DMAX1 ( DABS(DREAL(Z)), DABS(DIMAG(Z)) ) ! SI from "cdabs(z)" !!! onebyz = 1.0d0/z mz=zab + 1 mm=min(mz,12) p0=ap0(mm+1) p1=ap1(mm+1) q0=aq0(mm+1) q1=aq1(mm+1) z2 = onebyz*onebyz do 10 i=1,mm k=mm+1-i p0=p0*z2+ap0(k) p1=p1*z2+ap1(k) q0=q0*z2+aq0(k) 10 q1=q1*z2+aq1(k) a=cdsqrt(twobypi*onebyz) IF (DIMAG(Z).EQ.0.0D0) THEN IF (DREAL(Z).LT.0.0D0) A = -A ! correction by Sven Ivansson dec. 1990 END IF c Save four complex mults, by evaluating I* explicitly wor = q0*onebyz wor = dcmplx(-dimag(wor),dreal(wor)) h0 = a*(p0+wor) wor = z - pibyfour wor = dcmplx(-dimag(wor),dreal(wor)) h0 = h0*cdexp(wor) wor = q1*onebyz wor = dcmplx(-dimag(wor),dreal(wor)) h1 = a*(p1+wor) wor = z - thrpibyf wor = dcmplx(-dimag(wor),dreal(wor)) h1 = h1*cdexp(wor) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CASHANSPEC as cashan but with "common /inspec/" active C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine cashanspec(z,h0,h1) C C COMPUTES "inspec scaled" H0 AND H1 USING THE ASYMPTOTIC FORMULA C implicit real*8 (a-h,o-z) dimension ap0(13),ap1(13),aq0(13),aq1(13) complex*16 z,h0,h1,p0,p1,q0,q1,z2,a, onebyz, wor INCLUDE '../cpm/INSPEC.INC' parameter ( twobypi = .6366197723675814D+00 ) ! 2/PI parameter ( pibyfour = .7853981633974483D+00 ) ! PI/4 parameter ( thrpibyf = .2356194490192345D+01 ) ! 3*PI/4 DATA AP0 /.1D1,-.703125D-1,.112152099609375D0, ,-.5725014209747314D0,.6074042001273483D1, ,-.1100171402692467D3,.3038090510922384D4, ,-.1188384262567832D6,.6252951493434794D7, ,-.4259392165047667D9,.3646840080706554D11, ,-.3833534661393942D13,.4854014686852896D15/ DATA AP1 /.1D1,.1171875D0,-.144195556640625D0, ,.6765925884246826D0,-.6883914268109948D1, ,.1215978918765359D3,-.3302272294480852D4, ,.1276412726461746D6,-.6656367718817685D7, ,.4502786003050390D9,-.3833857520742787D11, ,.4011838599133194D13,-.5060568503314721D15/ DATA AQ0 /-.125D0,.732421875D-1,-.2271080017089844D0, ,.1727727502584457D1,-.2438052969955606D2, ,.5513358961220205D3,-.1825775547429317D5, ,.8328593040162890D6,-.5006958953198890D8, ,.3836255180230431D10,-.3649010818849831D12, ,.4218971570284093D14,-.5827244631566901D16/ DATA AQ1 /.375D0,-.1025390625D0,.2775764465332031D0, ,-.1993531733751297D1,.2724882731126854D2, ,-.6038440767050701D3,.1971837591223663D5, ,-.8902978767070675D6,.5310411010968520D8, ,-.4043620325107752D10,.3827011346598603D12, ,-.4406481417852275D14,.6065091351222692D16/ IF (DIMAG(Z).LT.0.0D0) THEN IF (DREAL(Z).LT.0.0D0) THEN WRITE(*,*) '%%% besbod.f, better accuracy by symmetry; Sven I stop' STOP 5760 END IF END IF zab = DMAX1 ( DABS(DREAL(Z)), DABS(DIMAG(Z)) ) ! SI from "cdabs(z)" !!! onebyz = 1.0d0/z mz=zab + 1 mm=min(mz,12) p0=ap0(mm+1) p1=ap1(mm+1) q0=aq0(mm+1) q1=aq1(mm+1) z2 = onebyz*onebyz do 10 i=1,mm k=mm+1-i p0=p0*z2+ap0(k) p1=p1*z2+ap1(k) q0=q0*z2+aq0(k) 10 q1=q1*z2+aq1(k) a=cdsqrt(twobypi*onebyz) IF (DIMAG(Z).EQ.0.0D0) THEN IF (DREAL(Z).LT.0.0D0) A = -A ! correction by Sven Ivansson dec. 1990 END IF c Save four complex mults, by evaluating I* explicitly wor = q0*onebyz wor = dcmplx(-dimag(wor),dreal(wor)) h0 = a*(p0+wor) wor = z - pibyfour wor = dcmplx(-dimag(wor),dreal(wor)) IF (IRKSPEC.EQ.0) THEN h0 = h0*cdexp(wor) ELSE h0 = h0*cdexp(wor+DLNSPEC) END IF wor = q1*onebyz wor = dcmplx(-dimag(wor),dreal(wor)) h1 = a*(p1+wor) wor = z - thrpibyf wor = dcmplx(-dimag(wor),dreal(wor)) IF (IRKSPEC.EQ.0) THEN h1 = h1*cdexp(wor) ELSE h1 = h1*cdexp(wor+DLNSPEC) END IF return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CSUMNEU C Computes Y0 and Y1 given J0 and J1, complex argument. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine csumneu(z,j0,j1,y0,y1) implicit real*8 (a-h,o-z) complex*16 z,j0,j1,y0,y1,z0,z2,f0,t0,f1,t1,a, onebyz parameter ( twobypi = .6366197723675814D+00 ) ! 2/PI parameter ( onebypi = .3183098861837907D+00 ) ! 1/PI onebyz = 1.0d0/z p0 = -.5772156649015329d0 p1 = p0 + 1.d0 z0 = 0.5d0*z a=cdlog(z0) y0=twobypi*(a*j0-p0) y1=(2*(a*j1-onebyz)-z0*(p0+p1))*onebypi z2=-z0*z0 f0=twobypi*z2 f1=0.5d0*onebypi*z0*z2 onebyi = 0.5d0 do 10 i=2,100 onebyi1 = 1.d0/(i+1) p0=p1 p1=p1 + onebyi t0=f0*p0 t1=f1*(p0+p1) y0=y0-t0 y1=y1-t1 IF ( DMAX1(DABS(DREAL(T0)),DABS(DIMAG(T0))) .LT. $1.D-13*DMAX1(DABS(DREAL(Y0)),DABS(DIMAG(Y0))) ) GOTO 20 C SI change from " if(cdabs(t0).lt.1.d-13*cdabs(y0)) goto 20 " !!! f0 = f0*z2*onebyi*onebyi f1 = f1*z2*onebyi*onebyi1 10 onebyi = onebyi1 stop 451 20 continue return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CSUMNEUSI C Computes "Y0 and Y1 except for the logarithmic terms and the C 1/z term", returned in y0reg and y1reg which are even entire C functions whose power-series have real coefficients, complex C argument. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine csumneusi(z,y0reg,y1reg) implicit real*8 (a-h,o-z) complex*16 z,y0reg,y1reg,z0,z2,f0,t0,f1,t1 parameter ( twobypi = .6366197723675814D+00 ) ! 2/PI parameter ( onebypi = .3183098861837907D+00 ) ! 1/PI p0 = -.5772156649015329d0 p1 = p0 + 1.d0 z0 = 0.5d0*z y0reg=twobypi*(-p0) y1reg=(-z0*(p0+p1))*onebypi z2=-z0*z0 f0=twobypi*z2 f1=0.5d0*onebypi*z0*z2 onebyi = 0.5d0 do 10 i=2,100 onebyi1 = 1.d0/(i+1) p0=p1 p1=p1 + onebyi t0=f0*p0 t1=f1*(p0+p1) y0reg=y0reg-t0 y1reg=y1reg-t1 IF ( DMAX1(DABS(DREAL(T0)),DABS(DIMAG(T0))) .LT. $1.D-13*DMAX1(DABS(DREAL(Y0REG)),DABS(DIMAG(Y0REG))) ) GOTO 20 C SI change from " if(cdabs(t0).lt.1.d-13*cdabs(y0reg)) goto 20 " !!! f0 = f0*z2*onebyi*onebyi f1 = f1*z2*onebyi*onebyi1 10 onebyi = onebyi1 stop 451 20 continue return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CDBES C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine cdbes(mmax1,za,b) C C COMPUTES THE BESSEL FUNTIONS OF ORDER 0,1,...,MMAX1-1 C AND COMPLEX ARGUMENT Z AND PUTS THEM IN THE ARRAY B(1),... C MMAX1<52; ABOUT 13 FIGURES ARE ACCURATE C USES SUBROUTINES CSUMBES CASBES C c 891209 IK: Argument ZA is now input parameter only. c Note that Z is changed below. implicit real*8 (a-h,o-z) complex*16 za,z,b(2),j0,j1,j2, twobyz parameter ( zimmax = 200.d0) ! Max ABS(Im(Z)) allowed c ! See comment first in file logical*2 ovfl, negaz z = za x = DMAX1 ( DABS(DREAL(Z)), DABS(DIMAG(Z)) ) ! SI from "cdabs(z)" !!! mx=x ovfl = abs(dimag(z)).gt.zimmax if(x.eq.0.d0.or.ovfl) then b(1) = dcmplx(1.d0,0.d0) if(ovfl) b(1) = dcmplx(0.d0,0.d0) do 1 i = 2, mmax1 1 b(i) = dcmplx(0.d0,0.d0) return end if twobyz = 2.0d0/z if(mx.gt.mmax1.and.mx.ge.20) then negaz = dreal(z).lt.0.d0 if(negaz) z = -z call casbes(z,j0,j1) if(negaz) z = -z if(negaz) j1 = -j1 b(1)=j0 c 891209 IK: Do not set B(2) when MMAX1 = 1 !! if(mmax1.ge.2) b(2)=j1 do 10 m = 3,mmax1 10 b(m) = (m-2)*twobyz*b(m-1) - b(m-2) return else mx=2*x-10 mm=max(mmax1,mx,3) call csumbes(mm-1,z,j2) ! J2 = Jmm-1(Z) call csumbes(mm-2,z,j1) ! J1 = Jmm-2(Z) if(mm.eq.mmax1) b(mm) = j2 if(mm-1.le.mmax1) b(mm-1) = j1 do 20 i=3,mm m=mm-i+1 j0 = twobyz*m*j1 - j2 j2 = j1 j1 = j0 20 if(m.le.mmax1) b(m) = j0 end if return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CSUMBES C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE CSUMBES(M,Z,JM) C C COMPUTES JM(Z) FOR COMPLEX Z USING THE SUM FORMULA C COMPLEX*16 Z,JM,TERM,Z2 JM=1.D0 TERM=1.D0 Z2=-.25D0*Z*Z do 10 i=1,100 kden = i*(i+m) term=term*z2 term=term/kden jm = jm + term 10 IF ( DMAX1(DABS(DREAL(TERM)),DABS(DIMAG(TERM))) .LT. $1.D-13 ) GOTO 20 C SI change from " 10 if(cdabs(term).lt.1.d-13) goto 20 " !!! 20 if(m.eq.0) return z2 = 0.5d0*z do 30 j=1,m jm = jm*z2 30 jm = dcmplx(dreal(jm)/j,dimag(jm)/j) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CASBES C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE CASBES(Z,J0,J1) C C COMPUTES J0 AND J1 FOR COMPLEX Z USING THE ASYMPTOTIC FORMULA C implicit real*8 (a-h,o-z) dimension ap0(13),ap1(13),aq0(13),aq1(13) complex*16 z,j0,j1,p0,p1,q0,q1,z2,y,s,c,a,onebyz parameter ( twobypi = .6366197723675814D+00 ) ! 2/PI parameter ( pibyfour = .7853981633974483D+00 ) ! PI/4 DATA AP0 /.1D1,-.703125D-1,.112152099609375D0, ,-.5725014209747314D0,.6074042001273483D1, ,-.1100171402692467D3,.3038090510922384D4, ,-.1188384262567832D6,.6252951493434794D7, ,-.4259392165047667D9,.3646840080706554D11, ,-.3833534661393942D13,.4854014686852896D15/ DATA AP1 /.1D1,.1171875D0,-.144195556640625D0, ,.6765925884246826D0,-.6883914268109948D1, ,.1215978918765359D3,-.3302272294480852D4, ,.1276412726461746D6,-.6656367718817685D7, ,.4502786003050390D9,-.3833857520742787D11, ,.4011838599133194D13,-.5060568503314721D15/ DATA AQ0 /-.125D0,.732421875D-1,-.2271080017089844D0, ,.1727727502584457D1,-.2438052969955606D2, ,.5513358961220205D3,-.1825775547429317D5, ,.8328593040162890D6,-.5006958953198890D8, ,.3836255180230431D10,-.3649010818849831D12, ,.4218971570284093D14,-.5827244631566901D16/ DATA AQ1 /.375D0,-.1025390625D0,.2775764465332031D0, ,-.1993531733751297D1,.2724882731126854D2, ,-.6038440767050701D3,.1971837591223663D5, ,-.8902978767070675D6,.5310411010968520D8, ,-.4043620325107752D10,.3827011346598603D12, ,-.4406481417852275D14,.6065091351222692D16/ DATA PI /3.141592653589793D0/ IF (DREAL(Z).LT.0.0D0) THEN WRITE(*,*) '%%% besbod.f, more accurate for "-z", J0 even; Sven I stop' STOP 5760 END IF P0=AP0(13) P1=AP1(13) Q0=AQ0(13) Q1=AQ1(13) onebyz = 1.0d0/z z2 = onebyz*onebyz DO 10 I=1,12 K=13-I P0=P0*Z2+AP0(K) P1=P1*Z2+AP1(K) Q0=Q0*Z2+AQ0(K) 10 Q1=Q1*Z2+AQ1(K) Y=Z-pibyfour S=CDSIN(Y) C=CDCOS(Y) a = cdsqrt(twobypi*onebyz) j0=a*(p0*c-q0*s*onebyz) j1=a*(p1*s+q1*c*onebyz) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C DHAN C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE DHAN(MMAX1,X,B,H) C C COMPUTES THE BESSEL AND HANKEL FUNTIONS OF ORDER 0,1,...,MMAX1-1 C AND REAL ARGUMENT X>0 AND PUTS THEM IN THE ARRAYS B(1),... AND H(1),... C MMAX1<52; ABOUT 12 FIGURES ARE ACCURATE C USES SUBROUTINES SUMBES, SUMNEU, AND ASBN C implicit real*8 (a-h,o-z) real*8 j0, j1, j2, b(2) complex*16 h(2) mx=x if(x.ne.0.d0) twobyx = 2.d0/x if(mx.gt.mmax1.and.mx.ge.20) then call asbn(x,j0,j1,y0,y1) b(1)=j0 if(mmax1.gt.1) b(2)=j1 do 10 m=3,mmax1 10 b(m)=twobyx*(m-2)*b(m-1) - b(m-2) else c 891209 IK: Avoid division by zero if(x.ne.0.d0) then c X is nonzero: c Use Taylor's formula and backward recursion w.r.t. order mx=2*x-10 mm=max(mmax1,mx,3) call sumbes(mm-1,x,j2) call sumbes(mm-2,x,j1) if(mm.le.mmax1) b(mm) = j2 if(mm-1.le.mmax1) b(mm-1) = j1 do 20 i=3,mm m=mm-i+1 j0 = twobyx*m*j1 - j2 j2 = j1 j1 = j0 20 if(m.le.mmax1) b(m) = j0 mx=x if(mx.lt.12) then c J0 and J1 are input data to SUMNEU call sumneu(x,j0,j1,y0,y1) else call asbn(x,j0,j1,y0,y1) endif else c X is zero c The bessel fcns of the second kind Y0, Y1,... and c the Hankel fcns have a pole at X = 0, and are thus c udefined. Their value is set to zero. b(1) = 1.d0 C "overflow" but choose symmetry in " J0 = 0.5d0 * ( H0(1) + H0(2) ) " H(1) = (1.0D0,0.0D0) !!! SI do 30 m = 3, mmax1 h(m) = dcmplx(0.d0,0.d0) 30 b(m) = 0.d0 return end if ! End X = 0 end if ! End use Taylor's formel and backw recursion c Here B(k) = Jk-1(X), k = 1,...,MMAX1 and c Y0 = Y0(X), Y1 = Y1(X) h(1) = dcmplx(b(1),y0) if(mmax1.ge.2) h(2) = dcmplx(b(2),y1) c When is the forward recursion below stable? do 40 m=3,mmax1 y2 = twobyx*(m-2)*y1 - y0 y0 = y1 y1 = y2 40 h(m) = dcmplx(b(m),y1) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUMBES C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE SUMBES(M,X,JM) C C COMPUTES JM(X) USING THE SUM FORMULA C implicit real*8 (a-h,o-z) real*8 jm JM=1.D0 TERM=1.D0 X2=-.25D0*X*X DO 10 I=1,100 TERM=TERM*X2/((M+I)*I) JM=JM+TERM 10 IF(DABS(TERM).LT.1.D-13) GOTO 20 20 IF(M.EQ.0) RETURN DO 30 J=1,M 30 JM=JM*X/(2*(M+1-J)) RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C ASBN C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE ASBN(X,J0,J1,Y0,Y1) C C COMPUTES J0,J1,Y0, AND Y1 USING THE ASYMPTOTIC FORMULA C implicit real*8 (a-h,o-z) real*8 j0, j1 dimension ap0(13),ap1(13),aq0(13),aq1(13) parameter ( twobypi = .6366197723675814D+00 ) ! 2/PI DATA AP0 /.1D1,-.703125D-1,.112152099609375D0, ,-.5725014209747314D0,.6074042001273483D1, ,-.1100171402692467D3,.3038090510922384D4, ,-.1188384262567832D6,.6252951493434794D7, ,-.4259392165047667D9,.3646840080706554D11, ,-.3833534661393942D13,.4854014686852896D15/ DATA AP1 /.1D1,.1171875D0,-.144195556640625D0, ,.6765925884246826D0,-.6883914268109948D1, ,.1215978918765359D3,-.3302272294480852D4, ,.1276412726461746D6,-.6656367718817685D7, ,.4502786003050390D9,-.3833857520742787D11, ,.4011838599133194D13,-.5060568503314721D15/ DATA AQ0 /-.125D0,.732421875D-1,-.2271080017089844D0, ,.1727727502584457D1,-.2438052969955606D2, ,.5513358961220205D3,-.1825775547429317D5, ,.8328593040162890D6,-.5006958953198890D8, ,.3836255180230431D10,-.3649010818849831D12, ,.4218971570284093D14,-.5827244631566901D16/ DATA AQ1 /.375D0,-.1025390625D0,.2775764465332031D0, ,-.1993531733751297D1,.2724882731126854D2, ,-.6038440767050701D3,.1971837591223663D5, ,-.8902978767070675D6,.5310411010968520D8, ,-.4043620325107752D10,.3827011346598603D12, ,-.4406481417852275D14,.6065091351222692D16/ DATA PI /3.141592653589793D0/ IF (X.LT.0.0D0) THEN WRITE(*,*) '%%% besbod.f, better accuracy by symmetry; Sven I stop' STOP 5760 END IF P0=AP0(13) P1=AP1(13) Q0=AQ0(13) Q1=AQ1(13) X2=1/(X*X) DO 10 I=1,12 K=13-I P0=P0*X2+AP0(K) P1=P1*X2+AP1(K) Q0=Q0*X2+AQ0(K) 10 Q1=Q1*X2+AQ1(K) Y=X-PI/4 S=DSIN(Y) C=DCOS(Y) a=dsqrt(twobypi/x) J0=A*(P0*C-Q0*S/X) J1=A*(P1*S+Q1*C/X) Y0=A*(P0*S+Q0*C/X) Y1=A*(-P1*C+Q1*S/X) RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUMNEU C Computes Y0 and Y1 given J0 and J1 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine sumneu(x,j0,j1,y0,y1) implicit real*8 (a-h,o-z) real*8 j0, j1 parameter ( twobypi = .6366197723675814D+00 ) ! 2/PI parameter ( onebypi = .3183098861837907D+00 ) ! 1/PI P0=-.5772156649015329D0 P1=P0+1 X0=X/2 A=DLOG(X0) Y0 = twobypi*(A*J0-P0) Y1 = onebypi*(2*A*J1-1/X0-X0*(P0+P1)) X2=-X0*X0 F0 = twobypi*x2 F1 = 0.5*onebypi*x0*x2 DO 10 I=2,100 P0=P1 P1=P1+1.D0/I T0=F0*P0 T1=F1*(P0+P1) Y0=Y0-T0 Y1=Y1-T1 IF(DABS(T0).LT.1.D-13*DABS(Y0)) GOTO 20 F0=F0*X2/(I*I) 10 F1=F1*X2/(I*(I+1)) STOP 451 20 CONTINUE RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C DBES C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE DBES(MMAX1,X,B) C C COMPUTES THE BESSEL FUNTIONS OF ORDER 0,1,...,MMAX1-1 C AND REAL ARGUMENT X AND PUTS THEM IN THE ARRAY B(1),... C MMAX1<52; ABOUT 13 FIGURES ARE ACCURATE C USES SUBROUTINES SUMBES AND ASBES C implicit real*8 (a-h,o-z) real*8 b(2),j0,j1,j2 X1=DABS(X) MX=X1 IF(MX.GT.MMAX1.AND.MX.GE.20) THEN CALL ASBES(X1,J0,J1) IF(X.LT.0.D0) J1=-J1 B(1)=J0 B(2)=J1 IF(MMAX1.LE.2) RETURN DO 10 M=3,MMAX1 10 B(M)=2*(M-2)*B(M-1)/X-B(M-2) ELSE c 891207, IK: Avoid division by zero. if(x.ne.0.d0) then c X is nonzero: WRITE(*,*) ' %%% ERROR IN DBES IN BESBOD HERE, NOT ', $ 'CORRECTED %%%' ! all components of b(1:mmax1) are not always updated STOP 100 ! Sven Ivansson 18/2 91, Klas Oestlund got error; USE BESRAT mx=2*x1-10 mm=max(mmax1,mx,3) call sumbes(mm-1,x,j2) call sumbes(mm-2,x,j1) if(mm.le.mmax1) b(mm) = j2 if(mm-1.le.mmax1) b(mm) = j1 twobyx = 2.d0/x do 20 i=3,mm m=mm-i+1 j0 = twobyx*m*j1 - j2 j2 = j1 j1 = j0 20 if(m.le.mmax1) b(m) = j0 else c X is zero: b(1) = 1.d0 do 30 m = 2,mmax1 30 b(m) = 0.d0 end if ! End X = 0 ENDIF RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C ASBES C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE ASBES(X,J0,J1) C C COMPUTES J0 AND J1 USING THE ASYMPTOTIC FORMULA C implicit real*8 (a-h,o-z) real*8 j0, j1 dimension ap0(13),ap1(13),aq0(13),aq1(13) parameter ( twobypi = .6366197723675814D+00 ) ! 2/PI parameter ( pibyfour = .7853981633974483D+00 ) ! PI/4 DATA AP0 /.1D1,-.703125D-1,.112152099609375D0, ,-.5725014209747314D0,.6074042001273483D1, ,-.1100171402692467D3,.3038090510922384D4, ,-.1188384262567832D6,.6252951493434794D7, ,-.4259392165047667D9,.3646840080706554D11, ,-.3833534661393942D13,.4854014686852896D15/ DATA AP1 /.1D1,.1171875D0,-.144195556640625D0, ,.6765925884246826D0,-.6883914268109948D1, ,.1215978918765359D3,-.3302272294480852D4, ,.1276412726461746D6,-.6656367718817685D7, ,.4502786003050390D9,-.3833857520742787D11, ,.4011838599133194D13,-.5060568503314721D15/ DATA AQ0 /-.125D0,.732421875D-1,-.2271080017089844D0, ,.1727727502584457D1,-.2438052969955606D2, ,.5513358961220205D3,-.1825775547429317D5, ,.8328593040162890D6,-.5006958953198890D8, ,.3836255180230431D10,-.3649010818849831D12, ,.4218971570284093D14,-.5827244631566901D16/ DATA AQ1 /.375D0,-.1025390625D0,.2775764465332031D0, ,-.1993531733751297D1,.2724882731126854D2, ,-.6038440767050701D3,.1971837591223663D5, ,-.8902978767070675D6,.5310411010968520D8, ,-.4043620325107752D10,.3827011346598603D12, ,-.4406481417852275D14,.6065091351222692D16/ IF (X.LT.0.0D0) THEN WRITE(*,*) '%%% besbod.f, better accuracy by symmetry; Sven I stop' STOP 5760 END IF P0=AP0(13) P1=AP1(13) Q0=AQ0(13) Q1=AQ1(13) X2=1/(X*X) DO 10 I=1,12 K=13-I P0=P0*X2+AP0(K) P1=P1*X2+AP1(K) Q0=Q0*X2+AQ0(K) 10 Q1=Q1*X2+AQ1(K) y=x-pibyfour S=DSIN(Y) C=DCOS(Y) a=dsqrt(twobypi/x) J0=A*(P0*C-Q0*S/X) J1=A*(P1*S+Q1*C/X) RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C DBES0BOD C Computes J0(Z) for real Z, by calling DBES C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ real*8 function dbes0bod(z) real*8 z,j0(2) C----------------------------------------------------------------- C About 13 figures are accurate C----------------------------------------------------------------- call dbes(1,z,j0) dbes0bod = j0(1) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CBES0BOD C Computes J0(Z) for complex Z, by calling CDBES C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ complex*16 function cbes0bod(z) complex*16 z,j0(2) C----------------------------------------------------------------- C About 13 figures are accurate C----------------------------------------------------------------- call cdbes(1,z,j0) cbes0bod = j0(1) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C DHAN0BOD C Computes H0(Z) for real Z, by calling DHAN C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ complex*16 function dhan0bod(z) real*8 z, j0(2) complex*16 h0(2) C----------------------------------------------------------------- C----------------------------------------------------------------- call dhan(1,z,j0,h0) dhan0bod = h0(1) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CHAN0BOD C Computes H0(Z) for complex Z, by calling CDHANSVEN (changed from cdhan) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ complex*16 function chan0bod(z) complex*16 z, h0(2) C----------------------------------------------------------------- C----------------------------------------------------------------- call cdhansven(1,z,h0) ! changed from "call cdhan(1,z,j0,h0)" chan0bod = h0(1) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CHAN0SPEC as chan0bod but with "common /inspec/" active C Computes "inspec scaled" H0(Z) for complex Z, by calling CDHANSPEC C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ complex*16 function chan0spec(z) complex*16 z, h0(2) C----------------------------------------------------------------- C----------------------------------------------------------------- call cdhanspec(1,z,h0) chan0spec = h0(1) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C DBES1BOD C Computes J1(Z) for real Z, by calling DBES C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ real*8 function dbes1bod(z) real*8 z,j1(2) C----------------------------------------------------------------- C About 13 figures are accurate C----------------------------------------------------------------- call dbes(2,z,j1) dbes1bod = j1(2) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CBES1BOD C Computes J1(Z) for complex Z, by calling CDBES C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ complex*16 function cbes1bod(z) complex*16 z,j1(2) C----------------------------------------------------------------- C About 13 figures are accurate C----------------------------------------------------------------- call cdbes(2,z,j1) cbes1bod = j1(2) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C DHAN1BOD C Computes H1(Z) for real Z, by calling DHAN C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ complex*16 function dhan1bod(z) real*8 z, j1(2) complex*16 h1(2) C----------------------------------------------------------------- C----------------------------------------------------------------- call dhan(2,z,j1,h1) dhan1bod = h1(2) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CHAN1BOD C Computes H1(Z) for complex Z, by calling CDHANSVEN (changed from cdhan) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ complex*16 function chan1bod(z) complex*16 z, h1(2) C----------------------------------------------------------------- C----------------------------------------------------------------- call cdhansven(2,z,h1) chan1bod = h1(2) return end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C CHAN1SPEC C Computes "inspec scaled" H1(Z) for complex Z, by calling CDHANSPEC C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ complex*16 function chan1spec(z) complex*16 z, h1(2) C----------------------------------------------------------------- C----------------------------------------------------------------- call cdhanspec(2,z,h1) chan1spec = h1(2) return end