c c c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c this is the end of the debugging code and the beginning of the c hankel function code proper. c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c c c subroutine hanks104(z,hanks,n,ifexpon) implicit real *8 (a-h,o-z) complex *16 z,hanks(*) c c This subroutine evaluates the first n+1 Hankel functions of the c argument z. The user also has the option of evaluating the c functions H_m(z) scaled by the (complex) coefficient e^{-i \cdot z}. c This option is provided via the parameter ifexpon (see below) c c c input parameters: c c z - the complex number for which the hankel functions c H_0, H_1 are to be evaluated c n - the highest order of any Hankel function to be evaluated c ifexpon - the integer parameter telling the subroutine whether c to calculate the actual values of the hankel functions, c or the values of Hankel functions scaled by e^{-i \cdot z}. c Permitted values: 0 and 1. c ifexpon = 1 will cause the subroutine to evaluate the Hankel functions c honestly c ifexpon = 0 will cause the subroutine to scale the Hankel functions c by e^{-i \cdot z}. c c output parameters: c c hanks - the first n+1 Hankel functions of the (complex) argument z. c Please note that hanks(1) is the Hankel function of order 0, c hanks(2) is the Hankel function of order 1, ..., hanks(n+1) c is the Hankel function of order n c c . . . evaluate the functions h0,h1 c call hank103(z,hanks(1),hanks(2),ifexpon) c c c conduct recursion c do 1200 i1=2,n c i=i1-1 c hanks(i1+1)=(2*i)/z*hanks(i1)-hanks(i1-1) 1200 continue c return end c c c c c subroutine hanks103(z,hanks,n,ifexpon) implicit real *8 (a-h,o-z) complex *16 z,hanks(*),cd,cdd c c This subroutine evaluates the first n+1 Hankel functions of the c argument z. The user also has the option of evaluating the c functions H_m(z) scaled by the (complex) coefficient e^{-i \cdot z}. c This option is provided via the parameter ifexpon (see below) c c c input parameters: c c z - the complex number for which the hankel functions c H_0, H_1 are to be evaluated c n - the highest order of any Hankel function to be evaluated c ifexpon - the integer parameter telling the subroutine whether c to calculate the actual values of the hankel functions, c or the values of Hankel functions scaled by e^{-i \cdot z}. c Permitted values: 0 and 1. c ifexpon = 1 will cause the subroutine to evaluate the Hankel functions c honestly c ifexpon = 0 will cause the subroutine to scale the Hankel functions c by e^{-i \cdot z}. c c output parameters: c c hanks - the first n+1 Hankel functions of the (complex) argument z. c Please note that hanks(1) is the Hankel function of order 0, c hanks(2) is the Hankel function of order 1, ..., hanks(n+1) c is the Hankel function of order n c c . . . evaluate the functions h0,h1 c call hank103(z,hanks(1),hanks(2),ifexpon) c c c conduct recursion c cd=2/z cdd=cd do 1200 i1=2,n c i=i1-1 c cccc hanks(i1+1)=(2*i)/z*hanks(i1)-hanks(i1-1) hanks(i1+1)=cdd*hanks(i1)-hanks(i1-1) c cdd=cdd+cd 1200 continue c return end c c c c c subroutine hank103(z,h0,h1,ifexpon) implicit real *8 (a-h,o-z) complex *16 z,h0,h1,h0u,h0r,h1u,h1r, 1 fj0,fj1,y0,y1,zu,zr,ima,ser2,ser3,z2, 2 cclog,cd data ima/(0.0d0,1.0d0)/,pi/0.31415926535897932D+01/ c c this subroutine evaluates the hankel functions H_0^1, H_1^1 c for an arbitrary user-specified complex number z. The user c also has the option of evaluating the functions h0, h1 c scaled by the (complex) coefficient e^{-i \cdot z}. This c subroutine is a modification of the subroutine hank102 c (see), different from the latter by having the parameter c ifexpon. Please note that the subroutine hank102 is in c turn a slightly accelerated version of the old hank101 c (see). The principal claim to fame of all three is that c they are valid on the whole complex plane, and are c reasonably accurate (14-digit relative accuracy) and c reasonably fast. Also, please note that all three have not c been carefully tested in the third quadrant (both x and y c negative); some sort of numerical trouble is possible c (though has not been observed) for LARGE z in the third c quadrant. c c input parameters: c c z - the complex number for which the hankel functions c H_0, H_1 are to be evaluated c ifexpon - the integer parameter telling the subroutine whether c to calculate the actual values of the hankel functions, c or the values of Hankel functions scaled by e^{-i \cdot z}. c Permitted values: 0 and 1. c ifexpon = 1 will cause the subroutine to evaluate the Hankel functions c honestly c ifexpon = 0 will cause the subroutine to scale the Hankel functions c by e^{-i \cdot z}. c c output parameters: c c h0, h1 - the said Hankel functions c c c . . . if z in the upper half-plane - act accordingly c if(imag(z) .lt. 0) goto 1400 call hank103u(z,ier,h0,h1,ifexpon) return 1400 continue c c if z is in the right lower quadrant - act accordingly c if(real(z) .lt. 0) goto 2000 call hank103r(z,ier,h0,h1,ifexpon) return 2000 continue c c z is in the left lower quadrant. compute c h0, h1 at the points zu, zr obtained from z by reflection c in the x and y axis, respectively c zu=dconjg(z) zr=-zu c call hank103u(zu,ier,h0u,h1u,ifexpon) call hank103r(zr,ier,h0r,h1r,ifexpon) if(ifexpon .eq. 1) goto 3000 subt=abs(imag(zu)) cd=exp(ima*zu-subt) h0u=h0u*cd h1u=h1u*cd cd=exp(ima*zr-subt) h0r=h0r*cd h1r=h1r*cd 3000 continue c c compute the functions j0, j1, y0, y1 c at the point zr c half=1 half=half/2 y0=(h0u+h0r)*half/ima fj0=-(h0u-h0r)*half c y1=-(h1u-h1r)*half/ima fj1=(h1u+h1r)*half c c finally, compute h0, h1 c c . . . calculate ser2, ser3 c z2=-dconjg(z) cclog=cdlog(z2) ser2=y0-fj0*2/pi*cclog ser3=y1-fj1*2/pi*cclog c c reflect all of these in the imaginary axis c fj0=dconjg(fj0) fj1=-dconjg(fj1) c ser2=dconjg(ser2) ser3=-dconjg(ser3) c c reconstitute y0, y1 c cclog=cdlog(z) y0=ser2+fj0*2/pi*cclog y1=ser3+fj1*2/pi*cclog c h0=fj0+ima*y0 h1=fj1+ima*y1 if(ifexpon .eq. 1) return cd=exp(-ima*z+subt) h0=h0*cd h1=h1*cd return end c c c c c subroutine hank103u(z,ier,h0,h1,ifexpon) implicit real *8 (a-h,o-z) complex *16 z,ima,cd,h0,h1,ccex,zzz9 real *8 c0p1(70) real *8 c1p1(70) real *8 c0p2(62) real *8 c1p2(62) data ima/(0.0d0,1.0d0)/ c c this subroutine evaluates the hankel functions H_0^1, H_1^1 c for a user-specified complex number z in the upper half-plane. c it is reasonably accurate (14-digit relative accuracy) c and reasonably fast. c c c input parameters: c c z - the complex number for which the hankel functions c H_0, H_1 are to be evaluated c c output parameters: c c ier - error return code. c ier=0 means successful conclusion c ier=4 means that z is not in the upper half-plane c h0, h1 - the said Hankel functions c data c0p1/ 1 -.6619836118357782D-12, -.6619836118612709D-12, 2 -.7307514264754200D-21, 0.3928160926261892D-10, 3 0.5712712520172854D-09, -.5712712519967086D-09, 4 -.1083820384008718D-07, -.1894529309455499D-18, 5 0.7528123700585197D-07, 0.7528123700841491D-07, 6 0.1356544045548053D-16, -.8147940452202855D-06, 7 -.3568198575016769D-05, 0.3568198574899888D-05, 8 0.2592083111345422D-04, 0.4209074870019400D-15, 9 -.7935843289157352D-04, -.7935843289415642D-04, a -.6848330800445365D-14, 0.4136028298630129D-03, 1 0.9210433149997867D-03, -.9210433149680665D-03, 2 -.3495306809056563D-02, -.6469844672213905D-13, 3 0.5573890502766937D-02, 0.5573890503000873D-02, 4 0.3767341857978150D-12, -.1439178509436339D-01, 5 -.1342403524448708D-01, 0.1342403524340215D-01, 6 0.8733016209933828D-02, 0.1400653553627576D-11, 7 0.2987361261932706D-01, 0.2987361261607835D-01, 8 -.3388096836339433D-11, -.1690673895793793D+00, 9 0.2838366762606121D+00, -.2838366762542546D+00, a 0.7045107746587499D+00, -.5363893133864181D-11, 1 -.7788044738211666D+00, -.7788044738130360D+00, 2 0.5524779104964783D-11, 0.1146003459721775D+01, 3 0.6930697486173089D+00, -.6930697486240221D+00, 4 -.7218270272305891D+00, 0.3633022466839301D-11, 5 0.3280924142354455D+00, 0.3280924142319602D+00, 6 -.1472323059106612D-11, -.2608421334424268D+00, 7 -.9031397649230536D-01, 0.9031397649339185D-01, 8 0.5401342784296321D-01, -.3464095071668884D-12, 9 -.1377057052946721D-01, -.1377057052927901D-01, a 0.4273263742980154D-13, 0.5877224130705015D-02, 1 0.1022508471962664D-02, -.1022508471978459D-02, 2 -.2789107903871137D-03, 0.2283984571396129D-14, 3 0.2799719727019427D-04, 0.2799719726970900D-04, 4 -.3371218242141487D-16, -.3682310515545645D-05, 5 -.1191412910090512D-06, 0.1191412910113518D-06/ c data c1p1/ 1 0.4428361927253983D-12, -.4428361927153559D-12, 2 -.2575693161635231D-10, -.2878656317479645D-21, 3 0.3658696304107867D-09, 0.3658696304188925D-09, 4 0.7463138750413651D-19, -.6748894854135266D-08, 5 -.4530098210372099D-07, 0.4530098210271137D-07, 6 0.4698787882823243D-06, 0.5343848349451927D-17, 7 -.1948662942158171D-05, -.1948662942204214D-05, 8 -.1658085463182409D-15, 0.1316906100496570D-04, 9 0.3645368564036497D-04, -.3645368563934748D-04, a -.1633458547818390D-03, -.2697770638600506D-14, 1 0.2816784976551660D-03, 0.2816784976676616D-03, 2 0.2548673351180060D-13, -.6106478245116582D-03, 3 0.2054057459296899D-03, -.2054057460218446D-03, 4 -.6254962367291260D-02, 0.1484073406594994D-12, 5 0.1952900562500057D-01, 0.1952900562457318D-01, 6 -.5517611343746895D-12, -.8528074392467523D-01, 7 -.1495138141086974D+00, 0.1495138141099772D+00, 8 0.4394907314508377D+00, -.1334677126491326D-11, 9 -.1113740586940341D+01, -.1113740586937837D+01, a 0.2113005088866033D-11, 0.1170212831401968D+01, 1 0.1262152242318805D+01, -.1262152242322008D+01, 2 -.1557810619605511D+01, 0.2176383208521897D-11, 3 0.8560741701626648D+00, 0.8560741701600203D+00, 4 -.1431161194996653D-11, -.8386735092525187D+00, 5 -.3651819176599290D+00, 0.3651819176613019D+00, 6 0.2811692367666517D+00, -.5799941348040361D-12, 7 -.9494630182937280D-01, -.9494630182894480D-01, 8 0.1364615527772751D-12, 0.5564896498129176D-01, 9 0.1395239688792536D-01, -.1395239688799950D-01, a -.5871314703753967D-02, 0.1683372473682212D-13, 1 0.1009157100083457D-02, 0.1009157100077235D-02, 2 -.8997331160162008D-15, -.2723724213360371D-03, 3 -.2708696587599713D-04, 0.2708696587618830D-04, 4 0.3533092798326666D-05, -.1328028586935163D-16, 5 -.1134616446885126D-06, -.1134616446876064D-06/ c data c0p2/ 1 0.5641895835516786D+00, -.5641895835516010D+00, 2 -.3902447089770041D-09, -.3334441074447365D-11, 3 -.7052368835911731D-01, -.7052368821797083D-01, 4 0.1957299315085370D-08, -.3126801711815631D-06, 5 -.3967331737107949D-01, 0.3967327747706934D-01, 6 0.6902866639752817D-04, 0.3178420816292497D-06, 7 0.4080457166061280D-01, 0.4080045784614144D-01, 8 -.2218731025620065D-04, 0.6518438331871517D-02, 9 0.9798339748600499D-01, -.9778028374972253D-01, a -.3151825524811773D+00, -.7995603166188139D-03, 1 0.1111323666639636D+01, 0.1116791178994330D+01, 2 0.1635711249533488D-01, -.8527067497983841D+01, 3 -.2595553689471247D+02, 0.2586942834408207D+02, 4 0.1345583522428299D+03, 0.2002017907999571D+00, 5 -.3086364384881525D+03, -.3094609382885628D+03, 6 -.1505974589617013D+01, 0.1250150715797207D+04, 7 0.2205210257679573D+04, -.2200328091885836D+04, 8 -.6724941072552172D+04, -.7018887749450317D+01, 9 0.8873498980910335D+04, 0.8891369384353965D+04, a 0.2008805099643591D+02, -.2030681426035686D+05, 1 -.2010017782384992D+05, 0.2006046282661137D+05, 2 0.3427941581102808D+05, 0.3432892927181724D+02, 3 -.2511417407338804D+05, -.2516567363193558D+05, 4 -.3318253740485142D+02, 0.3143940826027085D+05, 5 0.1658466564673543D+05, -.1654843151976437D+05, 6 -.1446345041326510D+05, -.1645433213663233D+02, 7 0.5094709396573681D+04, 0.5106816671258367D+04, 8 0.3470692471612145D+01, -.2797902324245621D+04, 9 -.5615581955514127D+03, 0.5601021281020627D+03, a 0.1463856702925587D+03, 0.1990076422327786D+00, 1 -.9334741618922085D+01, -.9361368967669095D+01/ c data c1p2/ 1 -.5641895835446003D+00, -.5641895835437973D+00, 2 0.3473016376419171D-10, -.3710264617214559D-09, 3 0.2115710836381847D+00, -.2115710851180242D+00, 4 0.3132928887334847D-06, 0.2064187785625558D-07, 5 -.6611954881267806D-01, -.6611997176900310D-01, 6 -.3386004893181560D-05, 0.7146557892862998D-04, 7 -.5728505088320786D-01, 0.5732906930408979D-01, 8 -.6884187195973806D-02, -.2383737409286457D-03, 9 0.1170452203794729D+00, 0.1192356405185651D+00, a 0.8652871239920498D-02, -.3366165876561572D+00, 1 -.1203989383538728D+01, 0.1144625888281483D+01, 2 0.9153684260534125D+01, 0.1781426600949249D+00, 3 -.2740411284066946D+02, -.2834461441294877D+02, 4 -.2192611071606340D+01, 0.1445470231392735D+03, 5 0.3361116314072906D+03, -.3270584743216529D+03, 6 -.1339254798224146D+04, -.1657618537130453D+02, 7 0.2327097844591252D+04, 0.2380960024514808D+04, 8 0.7760611776965994D+02, -.7162513471480693D+04, 9 -.9520608696419367D+04, 0.9322604506839242D+04, a 0.2144033447577134D+05, 0.2230232555182369D+03, 1 -.2087584364240919D+05, -.2131762020653283D+05, 2 -.3825699231499171D+03, 0.3582976792594737D+05, 3 0.2642632405857713D+05, -.2585137938787267D+05, 4 -.3251446505037506D+05, -.3710875194432116D+03, 5 0.1683805377643986D+05, 0.1724393921722052D+05, 6 0.1846128226280221D+03, -.1479735877145448D+05, 7 -.5258288893282565D+04, 0.5122237462705988D+04, 8 0.2831540486197358D+04, 0.3905972651440027D+02, 9 -.5562781548969544D+03, -.5726891190727206D+03, a -.2246192560136119D+01, 0.1465347141877978D+03, 1 0.9456733342595993D+01, -.9155767836700837D+01/ c c if the user-specified z is in the lower half-plane c - bomb out c ier=0 com=z if(imag(z) .ge. 0) goto 1200 ier=4 return 1200 continue c done=1 thresh1=1**2 thresh2=3.7**2 thresh3=20**2 c c check if if the user-specified z is in one of the c intermediate regimes c d=z*dconjg(z) if( (d .lt. thresh1) .or. (d .gt. thresh3) ) goto 3000 c c the user-specified z is in one of the intermediate regimes. c act accordingly c c if(d .gt. thresh2) goto 2000 c c z is in the first intermediate regime: its absolute value is c between 1 and 3.7. act accordingly c c . . . evaluate the expansion c cd=done/cdsqrt(z) c ccex=cd if(ifexpon .eq. 1) ccex=ccex*cdexp(ima*z) c zzz9=z**9 m=35 call hank103p(c0p1,m,cd,h0) h0=h0*ccex * zzz9 c call hank103p(c1p1,m,cd,h1) h1=h1*ccex * zzz9 return 2000 continue c c z is in the second intermediate regime: its absolute value is c between 3.7 and 20. act accordingly. c cd=done/cdsqrt(z) c ccex=cd if(ifexpon .eq. 1) ccex=ccex*cdexp(ima*z) m=31 call hank103p(c0p2,m,cd,h0) h0=h0*ccex c m=31 call hank103p(c1p2,m,cd,h1) h1=h1*ccex return 3000 continue c c z is either in the local regime or the asymptotic one. c if it is in the local regime - act accordingly. c if(d .gt. 50.d0) goto 4000 call hank103l(z,h0,h1,ifexpon) return c c z is in the asymptotic regime. act accordingly. c 4000 continue call hank103a(z,h0,h1,ifexpon) return end c c c c subroutine hank103p(p,m,z,f) implicit real *8 (a-h,o-z) complex *16 p(1),z,f c c evaluate a polynomial at a point c f=p(m) do 1200 i=m-1,1,-1 f=f*z+p(i) 1200 continue return end c c c c c subroutine hank103a(z,h0,h1,ifexpon) implicit real *8 (a-h,o-z) dimension p(18),q(18),p1(18),q1(18) complex *16 z,zinv,pp,qq,ima,h0,h1,pp1,qq1, 1 cccexp,cdd,cdumb,zinv22 data ima/(0.0d0,1.0d0)/,pi/0.31415926535897932D+01/, 1 done/1.0d0/,cdumb/ 2 (0.70710678118654757D+00,-.70710678118654746D+00)/ c data p/ 1 0.1000000000000000D+01, -.7031250000000000D-01, 2 0.1121520996093750D+00, -.5725014209747314D+00, 3 0.6074042001273483D+01, -.1100171402692467D+03, 4 0.3038090510922384D+04, -.1188384262567833D+06, 5 0.6252951493434797D+07, -.4259392165047669D+09, 6 0.3646840080706556D+11, -.3833534661393944D+13, 7 0.4854014686852901D+15, -.7286857349377657D+17, 8 0.1279721941975975D+20, -.2599382102726235D+22, 9 0.6046711487532401D+24, -.1597065525294211D+27/ c data q/ 1 -.1250000000000000D+00, 0.7324218750000000D-01, 2 -.2271080017089844D+00, 0.1727727502584457D+01, 3 -.2438052969955606D+02, 0.5513358961220206D+03, 4 -.1825775547429317D+05, 0.8328593040162893D+06, 5 -.5006958953198893D+08, 0.3836255180230434D+10, 6 -.3649010818849834D+12, 0.4218971570284096D+14, 7 -.5827244631566907D+16, 0.9476288099260110D+18, 8 -.1792162323051699D+21, 0.3900121292034000D+23, 9 -.9677028801069847D+25, 0.2715581773544907D+28/ data p1/ 1 0.1000000000000000D+01, 0.1171875000000000D+00, 2 -.1441955566406250D+00, 0.6765925884246826D+00, 3 -.6883914268109947D+01, 0.1215978918765359D+03, 4 -.3302272294480852D+04, 0.1276412726461746D+06, 5 -.6656367718817687D+07, 0.4502786003050393D+09, 6 -.3833857520742789D+11, 0.4011838599133198D+13, 7 -.5060568503314726D+15, 0.7572616461117957D+17, 8 -.1326257285320556D+20, 0.2687496750276277D+22, 9 -.6238670582374700D+24, 0.1644739123064188D+27/ c data q1/ 1 0.3750000000000000D+00, -.1025390625000000D+00, 2 0.2775764465332031D+00, -.1993531733751297D+01, 3 0.2724882731126854D+02, -.6038440767050702D+03, 4 0.1971837591223663D+05, -.8902978767070679D+06, 5 0.5310411010968522D+08, -.4043620325107754D+10, 6 0.3827011346598606D+12, -.4406481417852279D+14, 7 0.6065091351222699D+16, -.9833883876590680D+18, 8 0.1855045211579829D+21, -.4027994121281017D+23, 9 0.9974783533410457D+25, -.2794294288720121D+28/ c c evaluate the asymptotic expansion for h0,h1 at c the user-supplied point z, provided it is not c in the fourth quadrant c m=10 zinv=done/z c pp=p(m) pp1=p1(m) zinv22=zinv**2 c qq=q(m) qq1=q1(m) c do 1600 i=m-1,1,-1 pp=pp* zinv22+p(i) pp1=pp1* zinv22+p1(i) qq=qq* zinv22+q(i) qq1=qq1* zinv22+q1(i) 1600 continue c qq=qq*zinv qq1=qq1*zinv c cccexp=1 if(ifexpon .eq. 1) cccexp=cdexp(ima*z) c cdd=cdsqrt(2/pi*zinv) c h0=pp+ima*qq h0=cdd*cdumb*cccexp * h0 c h1=pp1+ima*qq1 h1=-cdd*cccexp*cdumb* h1*ima c return end c c c c c subroutine hank103l(z,h0,h1,ifexpon) implicit real *8 (a-h,o-z) dimension cj0(16),cj1(16),ser2(16),ser2der(16) complex *16 z,fj0,fj1,y0,y1,h0,h1,z2,cd,ima,cdddlog c data gamma/0.5772156649015328606d+00/ data ima/(0.0d0,1.0d0)/,pi/0.31415926535897932D+01/, 1 two/2.0d0/ c c this subroutine evaluates the hankel functions H_0^1, H_1^1 c for a user-specified complex number z in the local regime, c i. e. for cdabs(z) < 1 in the upper half-plane, c and for cdabs(z) < 4 in the lower half-plane, c it is reasonably accurate (14-digit relative accuracy) and c reasonably fast. c c input parameters: c c z - the complex number for which the hankel functions c H_0, H_1 are to be evaluated c c output parameters: c c h0, h1 - the said Hankel functions c data cj0/ 1 0.1000000000000000D+01, -.2500000000000000D+00, 2 0.1562500000000000D-01, -.4340277777777778D-03, 3 0.6781684027777778D-05, -.6781684027777778D-07, 4 0.4709502797067901D-09, -.2402807549524439D-11, 5 0.9385966990329841D-14, -.2896903392077112D-16, 6 0.7242258480192779D-19, -.1496334396734045D-21, 7 0.2597802772107717D-24, -.3842903509035085D-27, 8 0.4901662639075363D-30, -.5446291821194848D-33/ data cj1/ 1 -.5000000000000000D+00, 0.6250000000000000D-01, 2 -.2604166666666667D-02, 0.5425347222222222D-04, 3 -.6781684027777778D-06, 0.5651403356481481D-08, 4 -.3363930569334215D-10, 0.1501754718452775D-12, 5 -.5214426105738801D-15, 0.1448451696038556D-17, 6 -.3291935672814899D-20, 0.6234726653058522D-23, 7 -.9991549123491221D-26, 0.1372465538941102D-28, 8 -.1633887546358454D-31, 0.1701966194123390D-34/ data ser2/ 1 0.2500000000000000D+00, -.2343750000000000D-01, 2 0.7957175925925926D-03, -.1412850839120370D-04, 3 0.1548484519675926D-06, -.1153828185281636D-08, 4 0.6230136717695511D-11, -.2550971742728932D-13, 5 0.8195247730999099D-16, -.2121234517551702D-18, 6 0.4518746345057852D-21, -.8061529302289970D-24, 7 0.1222094716680443D-26, -.1593806157473552D-29, 8 0.1807204342667468D-32, -.1798089518115172D-35/ data ser2der/ 1 0.5000000000000000D+00, -.9375000000000000D-01, 2 0.4774305555555556D-02, -.1130280671296296D-03, 3 0.1548484519675926D-05, -.1384593822337963D-07, 4 0.8722191404773715D-10, -.4081554788366291D-12, 5 0.1475144591579838D-14, -.4242469035103405D-17, 6 0.9941241959127275D-20, -.1934767032549593D-22, 7 0.3177446263369152D-25, -.4462657240925946D-28, 8 0.5421613028002404D-31, -.5753886457968550D-34/ c c evaluate j0, j1 c m=16 fj0=0 fj1=0 y0=0 y1=0 z2=z**2 cd=1 c do 1800 i=1,m fj0=fj0+cj0(i)*cd fj1=fj1+cj1(i)*cd y1=y1+ser2der(i)*cd cd=cd*z2 y0=y0+ser2(i)*cd 1800 continue fj1=-fj1*z c cdddlog=cdlog(z/two)+gamma y0=cdddlog*fj0+y0 y0=two/pi*y0 c y1=y1*z c y1=-cdddlog*fj1+fj0/z+y1 y1=-y1*two/pi c h0=fj0+ima*y0 h1=fj1+ima*y1 c if(ifexpon .eq. 1) return c cd=exp(-ima*z) h0=h0*cd h1=h1*cd c return end c c c c c subroutine hank103r(z,ier,h0,h1,ifexpon) implicit real *8 (a-h,o-z) complex *16 z,ima,cd,h0,h1,cccexp,cdd,zz18 real *8 c0p1(70) real *8 c1p1(70) real *8 c0p2(54) real *8 c1p2(62) data ima/(0.0d0,1.0d0)/ c c this subroutine evaluates the hankel functions H_0^1, H_1^1 c for a user-specified complex number z in the right lower c quadrant. it is reasonably accurate (14-digit relative c accuracy) and reasonably fast. c c c input parameters: c c z - the complex number for which the hankel functions c H_0, H_1 are to be evaluated c c output parameters: c c ier - error return code. c ier=0 means successful conclusion c ier=4 means that z is not in the right lower quadrant c h0, h1 - the said Hankel functions c data c0p1/ 1 -.4268441995428495D-23, 0.4374027848105921D-23, 2 0.9876152216238049D-23, -.1065264808278614D-20, 3 0.6240598085551175D-19, 0.6658529985490110D-19, 4 -.5107210870050163D-17, -.2931746613593983D-18, 5 0.1611018217758854D-15, -.1359809022054077D-15, 6 -.7718746693707326D-15, 0.6759496139812828D-14, 7 -.1067620915195442D-12, -.1434699000145826D-12, 8 0.3868453040754264D-11, 0.7061853392585180D-12, 9 -.6220133527871203D-10, 0.3957226744337817D-10, a 0.3080863675628417D-09, -.1154618431281900D-08, 1 0.7793319486868695D-08, 0.1502570745460228D-07, 2 -.1978090852638430D-06, -.7396691873499030D-07, 3 0.2175857247417038D-05, -.8473534855334919D-06, 4 -.1053381327609720D-04, 0.2042555121261223D-04, 5 -.4812568848956982D-04, -.1961519090873697D-03, 6 0.1291714391689374D-02, 0.9234422384950050D-03, 7 -.1113890671502769D-01, 0.9053687375483149D-03, 8 0.5030666896877862D-01, -.4923119348218356D-01, 9 0.5202355973926321D+00, -.1705244841954454D+00, a -.1134990486611273D+01, -.1747542851820576D+01, 1 0.8308174484970718D+01, 0.2952358687641577D+01, 2 -.3286074510100263D+02, 0.1126542966971545D+02, 3 0.6576015458463394D+02, -.1006116996293757D+03, 4 0.3216834899377392D+02, 0.3614005342307463D+03, 5 -.6653878500833375D+03, -.6883582242804924D+03, 6 0.2193362007156572D+04, 0.2423724600546293D+03, 7 -.3665925878308203D+04, 0.2474933189642588D+04, 8 0.1987663383445796D+04, -.7382586600895061D+04, 9 0.4991253411017503D+04, 0.1008505017740918D+05, a -.1285284928905621D+05, -.5153674821668470D+04, 1 0.1301656757246985D+05, -.4821250366504323D+04, 2 -.4982112643422311D+04, 0.9694070195648748D+04, 3 -.1685723189234701D+04, -.6065143678129265D+04, 4 0.2029510635584355D+04, 0.1244402339119502D+04, 5 -.4336682903961364D+03, 0.8923209875101459D+02/ c data c1p1/ 1 -.4019450270734195D-23, -.4819240943285824D-23, 2 0.1087220822839791D-20, 0.1219058342725899D-21, 3 -.7458149572694168D-19, 0.5677825613414602D-19, 4 0.8351815799518541D-18, -.5188585543982425D-17, 5 0.1221075065755962D-15, 0.1789261470637227D-15, 6 -.6829972121890858D-14, -.1497462301804588D-14, 7 0.1579028042950957D-12, -.9414960303758800D-13, 8 -.1127570848999746D-11, 0.3883137940932639D-11, 9 -.3397569083776586D-10, -.6779059427459179D-10, a 0.1149529442506273D-08, 0.4363087909873751D-09, 1 -.1620182360840298D-07, 0.6404695607668289D-08, 2 0.9651461037419628D-07, -.1948572160668177D-06, 3 0.6397881896749446D-06, 0.2318661930507743D-05, 4 -.1983192412396578D-04, -.1294811208715315D-04, 5 0.2062663873080766D-03, -.2867633324735777D-04, 6 -.1084309075952914D-02, 0.1227880935969686D-02, 7 0.2538406015667726D-03, -.1153316815955356D-01, 8 0.4520140008266983D-01, 0.5693944718258218D-01, 9 -.9640790976658534D+00, -.6517135574036008D+00, a 0.2051491829570049D+01, -.1124151010077572D+01, 1 -.3977380460328048D+01, 0.8200665483661009D+01, 2 -.7950131652215817D+01, -.3503037697046647D+02, 3 0.9607320812492044D+02, 0.7894079689858070D+02, 4 -.3749002890488298D+03, -.8153831134140778D+01, 5 0.7824282518763973D+03, -.6035276543352174D+03, 6 -.5004685759675768D+03, 0.2219009060854551D+04, 7 -.2111301101664672D+04, -.4035632271617418D+04, 8 0.7319737262526823D+04, 0.2878734389521922D+04, 9 -.1087404934318719D+05, 0.3945740567322783D+04, a 0.6727823761148537D+04, -.1253555346597302D+05, 1 0.3440468371829973D+04, 0.1383240926370073D+05, 2 -.9324927373036743D+04, -.6181580304530313D+04, 3 0.6376198146666679D+04, -.1033615527971958D+04, 4 -.1497604891055181D+04, 0.1929025541588262D+04, 5 -.4219760183545219D+02, -.4521162915353207D+03/ c data c0p2/ 1 0.5641895835569398D+00, -.5641895835321127D+00, 2 -.7052370223565544D-01, -.7052369923405479D-01, 3 -.3966909368581382D-01, 0.3966934297088857D-01, 4 0.4130698137268744D-01, 0.4136196771522681D-01, 5 0.6240742346896508D-01, -.6553556513852438D-01, 6 -.3258849904760676D-01, -.7998036854222177D-01, 7 -.3988006311955270D+01, 0.1327373751674479D+01, 8 0.6121789346915312D+02, -.9251865216627577D+02, 9 0.4247064992018806D+03, 0.2692553333489150D+04, a -.4374691601489926D+05, -.3625248208112831D+05, 1 0.1010975818048476D+07, -.2859360062580096D+05, 2 -.1138970241206912D+08, 0.1051097979526042D+08, 3 0.2284038899211195D+08, -.2038012515235694D+09, 4 0.1325194353842857D+10, 0.1937443530361381D+10, 5 -.2245999018652171D+11, -.5998903865344352D+10, 6 0.1793237054876609D+12, -.8625159882306147D+11, 7 -.5887763042735203D+12, 0.1345331284205280D+13, 8 -.2743432269370813D+13, -.8894942160272255D+13, 9 0.4276463113794564D+14, 0.2665019886647781D+14, a -.2280727423955498D+15, 0.3686908790553973D+14, 1 0.5639846318168615D+15, -.6841529051615703D+15, 2 0.9901426799966038D+14, 0.2798406605978152D+16, 3 -.4910062244008171D+16, -.5126937967581805D+16, 4 0.1387292951936756D+17, 0.1043295727224325D+16, 5 -.1565204120687265D+17, 0.1215262806973577D+17, 6 0.3133802397107054D+16, -.1801394550807078D+17, 7 0.4427598668012807D+16, 0.6923499968336864D+16/ c c data c1p2/ 1 -.5641895835431980D+00, -.5641895835508094D+00, 2 0.2115710934750869D+00, -.2115710923186134D+00, 3 -.6611607335011594D-01, -.6611615414079688D-01, 4 -.5783289433408652D-01, 0.5785737744023628D-01, 5 0.8018419623822896D-01, 0.8189816020440689D-01, 6 0.1821045296781145D+00, -.2179738973008740D+00, 7 0.5544705668143094D+00, 0.2224466316444440D+01, 8 -.8563271248520645D+02, -.4394325758429441D+02, 9 0.2720627547071340D+04, -.6705390850875292D+03, a -.3936221960600770D+05, 0.5791730432605451D+05, 1 -.1976787738827811D+06, -.1502498631245144D+07, 2 0.2155317823990686D+08, 0.1870953796705298D+08, 3 -.4703995711098311D+09, 0.3716595906453190D+07, 4 0.5080557859012385D+10, -.4534199223888966D+10, 5 -.1064438211647413D+11, 0.8612243893745942D+11, 6 -.5466017687785078D+12, -.8070950386640701D+12, 7 0.9337074941225827D+13, 0.2458379240643264D+13, 8 -.7548692171244579D+14, 0.3751093169954336D+14, 9 0.2460677431350039D+15, -.5991919372881911D+15, a 0.1425679408434606D+16, 0.4132221939781502D+16, 1 -.2247506469468969D+17, -.1269771078165026D+17, 2 0.1297336292749026D+18, -.2802626909791308D+17, 3 -.3467137222813017D+18, 0.4773955215582192D+18, 4 -.2347165776580206D+18, -.2233638097535785D+19, 5 0.5382350866778548D+19, 0.4820328886922998D+19, 6 -.1928978948099345D+20, 0.1575498747750907D+18, 7 0.3049162180215152D+20, -.2837046201123502D+20, 8 -.5429391644354291D+19, 0.6974653380104308D+20, 9 -.5322120857794536D+20, -.6739879079691706D+20, a 0.6780343087166473D+20, 0.1053455984204666D+20, 1 -.2218784058435737D+20, 0.1505391868530062D+20/ c c if z is not in the right lower quadrant - bomb out c ier=0 if( (real(z) .ge. 0) .and. (imag(z) .le. 0) ) goto 1400 ier=4 return 1400 continue c done=1 thresh1=4**2 thresh2=8**2 thresh3=20**2 c c check if if the user-specified z is in one of the c intermediate regimes c d=z*dconjg(z) if( (d .lt. thresh1) .or. (d .gt. thresh3) ) goto 3000 c c if the user-specified z is in the first intermediate regime c (i.e. if its absolute value is between 4 and 8), act accordingly c if(d .gt. thresh2) goto 2000 c cccexp=1 if(ifexpon .eq. 1) cccexp=cdexp(ima*z) cdd=done/cdsqrt(z) cd=done/z zz18=z**18 m=35 call hank103p(c0p1,m,cd,h0) h0=h0*cdd*cccexp*zz18 c call hank103p(c1p1,m,cd,h1) h1=h1*cdd*cccexp*zz18 return 2000 continue c c z is in the second intermediate regime (i.e. its c absolute value is between 8 and 20). act accordingly. c cd=done/z cdd=sqrt(cd) cccexp=1 if(ifexpon .eq. 1) cccexp=cdexp(ima*z) m=27 c call hank103p(c0p2,m,cd,h0) h0=h0*cccexp*cdd c m=31 call hank103p(c1p2,m,cd,h1) h1=h1*cccexp*cdd return 3000 continue c c c z is either in the local regime or the asymptotic one. c if it is in the local regime - act accordingly. c if(d .gt. 50.d0) goto 4000 call hank103l(z,h0,h1,ifexpon) return c c z is in the asymptotic regime. act accordingly. c 4000 continue call hank103a(z,h0,h1,ifexpon) return end