d9b0mp.f 14 KB
subroutine d9b0mp (x, ampl, theta)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
c
c evaluate the modulus and phase for the bessel j0 and y0 functions.
c
      double precision x, ampl, theta, bm0cs(37), bt02cs(39),
     1  bm02cs(40), bth0cs(44), xmax, pi4, z, d1mach, dcsevl,
     2  dsqrt
      external d1mach, dcsevl, initds
c
c series for bm0        on the interval  1.56250e-02 to  6.25000e-02
c                                        with weighted error   4.40e-32
c                                         log weighted error  31.36
c                               significant figures required  30.02
c                                    decimal places required  32.14
c
      data bm0 cs(  1) / +.9211656246 8277427125 7376773018 2 d-1      /
      data bm0 cs(  2) / -.1050590997 2719051024 8071637175 5 d-2      /
      data bm0 cs(  3) / +.1470159840 7687597540 5639285095 2 d-4      /
      data bm0 cs(  4) / -.5058557606 0385542233 4792932770 2 d-6      /
      data bm0 cs(  5) / +.2787254538 6324441766 3035613788 1 d-7      /
      data bm0 cs(  6) / -.2062363611 7809148026 1884101897 3 d-8      /
      data bm0 cs(  7) / +.1870214313 1388796751 3817259626 1 d-9      /
      data bm0 cs(  8) / -.1969330971 1356362002 4173077782 5 d-10     /
      data bm0 cs(  9) / +.2325973793 9992754440 1250881805 2 d-11     /
      data bm0 cs( 10) / -.3009520344 9382502728 5122473448 2 d-12     /
      data bm0 cs( 11) / +.4194521333 8506691814 7120676864 6 d-13     /
      data bm0 cs( 12) / -.6219449312 1884458259 7326742956 4 d-14     /
      data bm0 cs( 13) / +.9718260411 3360684696 0176588526 9 d-15     /
      data bm0 cs( 14) / -.1588478585 7010752073 6663596693 7 d-15     /
      data bm0 cs( 15) / +.2700072193 6713088900 8621732445 8 d-16     /
      data bm0 cs( 16) / -.4750092365 2340089924 7750478677 3 d-17     /
      data bm0 cs( 17) / +.8615128162 6043708731 9170374656 0 d-18     /
      data bm0 cs( 18) / -.1605608686 9561448157 4560270335 9 d-18     /
      data bm0 cs( 19) / +.3066513987 3144829751 8853980159 9 d-19     /
      data bm0 cs( 20) / -.5987764223 1939564306 9650561706 6 d-20     /
      data bm0 cs( 21) / +.1192971253 7482483064 8906984106 6 d-20     /
      data bm0 cs( 22) / -.2420969142 0448054894 8468258133 3 d-21     /
      data bm0 cs( 23) / +.4996751760 5106164533 7100287999 9 d-22     /
      data bm0 cs( 24) / -.1047493639 3511585100 9504051199 9 d-22     /
      data bm0 cs( 25) / +.2227786843 7974681010 4818346666 6 d-23     /
      data bm0 cs( 26) / -.4801813239 3981628623 7054293333 3 d-24     /
      data bm0 cs( 27) / +.1047962723 4709599564 7699626666 6 d-24     /
      data bm0 cs( 28) / -.2313858165 6786153251 0126080000 0 d-25     /
      data bm0 cs( 29) / +.5164823088 4626742116 3519999999 9 d-26     /
      data bm0 cs( 30) / -.1164691191 8500653895 2540159999 9 d-26     /
      data bm0 cs( 31) / +.2651788486 0433192829 5833600000 0 d-27     /
      data bm0 cs( 32) / -.6092559503 8257284976 9130666666 6 d-28     /
      data bm0 cs( 33) / +.1411804686 1442593080 3882666666 6 d-28     /
      data bm0 cs( 34) / -.3298094961 2317372457 5061333333 3 d-29     /
      data bm0 cs( 35) / +.7763931143 0740650317 1413333333 3 d-30     /
      data bm0 cs( 36) / -.1841031343 6614584784 2133333333 3 d-30     /
      data bm0 cs( 37) / +.4395880138 5943107371 0079999999 9 d-31     /
c
c series for bth0       on the interval  0.          to  1.56250e-02
c                                        with weighted error   2.66e-32
c                                         log weighted error  31.57
c                               significant figures required  30.67
c                                    decimal places required  32.40
c
      data bth0cs(  1) / -.2490178086 2128936717 7097937899 67 d+0     /
      data bth0cs(  2) / +.4855029960 9623749241 0486155354 85 d-3     /
      data bth0cs(  3) / -.5451183734 5017204950 6562735635 05 d-5     /
      data bth0cs(  4) / +.1355867305 9405964054 3774459299 03 d-6     /
      data bth0cs(  5) / -.5569139890 2227626227 5832184149 20 d-8     /
      data bth0cs(  6) / +.3260903182 4994335304 0042057194 68 d-9     /
      data bth0cs(  7) / -.2491880786 2461341125 2379038779 93 d-10    /
      data bth0cs(  8) / +.2344937742 0882520554 3524135648 91 d-11    /
      data bth0cs(  9) / -.2609653444 4310387762 1775747661 36 d-12    /
      data bth0cs( 10) / +.3335314042 0097395105 8699550149 23 d-13    /
      data bth0cs( 11) / -.4789000044 0572684646 7507705574 09 d-14    /
      data bth0cs( 12) / +.7595617843 6192215972 6425685452 48 d-15    /
      data bth0cs( 13) / -.1313155601 6891440382 7733974876 33 d-15    /
      data bth0cs( 14) / +.2448361834 5240857495 4268207383 55 d-16    /
      data bth0cs( 15) / -.4880572981 0618777683 2567619183 31 d-17    /
      data bth0cs( 16) / +.1032728502 9786316149 2237563612 04 d-17    /
      data bth0cs( 17) / -.2305763381 5057217157 0047445270 25 d-18    /
      data bth0cs( 18) / +.5404444300 1892693993 0171084837 65 d-19    /
      data bth0cs( 19) / -.1324069519 4366572724 1550328823 85 d-19    /
      data bth0cs( 20) / +.3378079562 1371970203 4247921247 22 d-20    /
      data bth0cs( 21) / -.8945762915 7111779003 0269262922 99 d-21    /
      data bth0cs( 22) / +.2451990688 9219317090 8999086514 05 d-21    /
      data bth0cs( 23) / -.6938842287 6866318680 1399331576 57 d-22    /
      data bth0cs( 24) / +.2022827871 4890138392 9463033377 91 d-22    /
      data bth0cs( 25) / -.6062850000 2335483105 7941953717 64 d-23    /
      data bth0cs( 26) / +.1864974896 4037635381 8237883962 70 d-23    /
      data bth0cs( 27) / -.5878373238 4849894560 2450365308 67 d-24    /
      data bth0cs( 28) / +.1895859144 7999563485 5311795035 13 d-24    /
      data bth0cs( 29) / -.6248197937 2258858959 2916207285 65 d-25    /
      data bth0cs( 30) / +.2101790168 4551024686 6386335290 74 d-25    /
      data bth0cs( 31) / -.7208430093 5209253690 8139339924 46 d-26    /
      data bth0cs( 32) / +.2518136389 2474240867 1564059767 46 d-26    /
      data bth0cs( 33) / -.8951804225 8785778806 1439459536 43 d-27    /
      data bth0cs( 34) / +.3235723747 9762298533 2562358685 87 d-27    /
      data bth0cs( 35) / -.1188301051 9855353657 0471441137 96 d-27    /
      data bth0cs( 36) / +.4430628690 7358104820 5792319417 31 d-28    /
      data bth0cs( 37) / -.1676100964 8834829495 7920101356 81 d-28    /
      data bth0cs( 38) / +.6429294692 1207466972 5323939660 88 d-29    /
      data bth0cs( 39) / -.2499226116 6978652421 2072136827 63 d-29    /
      data bth0cs( 40) / +.9839979429 9521955672 8282603553 18 d-30    /
      data bth0cs( 41) / -.3922037524 2408016397 9891316261 58 d-30    /
      data bth0cs( 42) / +.1581810703 0056522138 5906188456 92 d-30    /
      data bth0cs( 43) / -.6452550614 4890715944 3440983654 26 d-31    /
      data bth0cs( 44) / +.2661111136 9199356137 1770183463 67 d-31    /
c
c series for bm02       on the interval  0.          to  1.56250e-02
c                                        with weighted error   4.72e-32
c                                         log weighted error  31.33
c                               significant figures required  30.00
c                                    decimal places required  32.13
c
      data bm02cs(  1) / +.9500415145 2283813693 3086133556 0 d-1      /
      data bm02cs(  2) / -.3801864682 3656709917 4808156685 1 d-3      /
      data bm02cs(  3) / +.2258339301 0314811929 5182992722 4 d-5      /
      data bm02cs(  4) / -.3895725802 3722287647 3062141260 5 d-7      /
      data bm02cs(  5) / +.1246886416 5120816979 3099052972 5 d-8      /
      data bm02cs(  6) / -.6065949022 1025037798 0383505838 7 d-10     /
      data bm02cs(  7) / +.4008461651 4217469910 1527597104 5 d-11     /
      data bm02cs(  8) / -.3350998183 3980942184 6729879457 4 d-12     /
      data bm02cs(  9) / +.3377119716 5174173670 6326434199 6 d-13     /
      data bm02cs( 10) / -.3964585901 6350127005 6935629582 3 d-14     /
      data bm02cs( 11) / +.5286111503 8838572173 8793974473 5 d-15     /
      data bm02cs( 12) / -.7852519083 4508523136 5464024349 3 d-16     /
      data bm02cs( 13) / +.1280300573 3866822010 1163407344 9 d-16     /
      data bm02cs( 14) / -.2263996296 3914297762 8709924488 4 d-17     /
      data bm02cs( 15) / +.4300496929 6567903886 4641029047 7 d-18     /
      data bm02cs( 16) / -.8705749805 1325870797 4753545145 5 d-19     /
      data bm02cs( 17) / +.1865862713 9620951411 8144277205 0 d-19     /
      data bm02cs( 18) / -.4210482486 0930654573 4508697230 1 d-20     /
      data bm02cs( 19) / +.9956676964 2284009915 8162741784 2 d-21     /
      data bm02cs( 20) / -.2457357442 8053133596 0592147854 7 d-21     /
      data bm02cs( 21) / +.6307692160 7620315680 8735370705 9 d-22     /
      data bm02cs( 22) / -.1678773691 4407401426 9333117238 8 d-22     /
      data bm02cs( 23) / +.4620259064 6739044337 7087813608 7 d-23     /
      data bm02cs( 24) / -.1311782266 8603087322 3769340249 6 d-23     /
      data bm02cs( 25) / +.3834087564 1163028277 4792244027 6 d-24     /
      data bm02cs( 26) / -.1151459324 0777412710 7261329357 6 d-24     /
      data bm02cs( 27) / +.3547210007 5233385230 7697134521 3 d-25     /
      data bm02cs( 28) / -.1119218385 8150046462 6435594217 6 d-25     /
      data bm02cs( 29) / +.3611879427 6298378316 9840499425 7 d-26     /
      data bm02cs( 30) / -.1190687765 9133331500 9264176246 3 d-26     /
      data bm02cs( 31) / +.4005094059 4039681318 0247644953 6 d-27     /
      data bm02cs( 32) / -.1373169422 4522123905 9519391601 7 d-27     /
      data bm02cs( 33) / +.4794199088 7425315859 9649152643 7 d-28     /
      data bm02cs( 34) / -.1702965627 6241095840 0699447645 2 d-28     /
      data bm02cs( 35) / +.6149512428 9363300715 0357516132 4 d-29     /
      data bm02cs( 36) / -.2255766896 5818283499 4430023724 2 d-29     /
      data bm02cs( 37) / +.8399707509 2942994860 6165835320 0 d-30     /
      data bm02cs( 38) / -.3172997595 5626023555 6742393615 2 d-30     /
      data bm02cs( 39) / +.1215205298 8812985545 8333302651 4 d-30     /
      data bm02cs( 40) / -.4715852749 7544386930 1321056804 5 d-31     /
c
c series for bt02       on the interval  1.56250e-02 to  6.25000e-02
c                                        with weighted error   2.99e-32
c                                         log weighted error  31.52
c                               significant figures required  30.61
c                                    decimal places required  32.32
c
      data bt02cs(  1) / -.2454829521 3424597462 0504672493 24 d+0     /
      data bt02cs(  2) / +.1254412103 9084615780 7853317782 99 d-2     /
      data bt02cs(  3) / -.3125395041 4871522854 9734467095 71 d-4     /
      data bt02cs(  4) / +.1470977824 9940831164 4534269693 14 d-5     /
      data bt02cs(  5) / -.9954348893 7950033643 4688503511 58 d-7     /
      data bt02cs(  6) / +.8549316673 3203041247 5787113977 51 d-8     /
      data bt02cs(  7) / -.8698975952 6554334557 9855121791 92 d-9     /
      data bt02cs(  8) / +.1005209953 3559791084 5401010821 53 d-9     /
      data bt02cs(  9) / -.1282823060 1708892903 4836236855 44 d-10    /
      data bt02cs( 10) / +.1773170078 1805131705 6557504510 23 d-11    /
      data bt02cs( 11) / -.2617457456 9485577488 6362841809 25 d-12    /
      data bt02cs( 12) / +.4082835138 9972059621 9664812211 03 d-13    /
      data bt02cs( 13) / -.6675166823 9742720054 6067495542 61 d-14    /
      data bt02cs( 14) / +.1136576139 3071629448 3924695499 51 d-14    /
      data bt02cs( 15) / -.2005118962 0647160250 5592664121 17 d-15    /
      data bt02cs( 16) / +.3649797879 4766269635 7205914641 06 d-16    /
      data bt02cs( 17) / -.6830963756 4582303169 3558437888 00 d-17    /
      data bt02cs( 18) / +.1310758314 5670756620 0571042679 46 d-17    /
      data bt02cs( 19) / -.2572336310 1850607778 7571306495 99 d-18    /
      data bt02cs( 20) / +.5152165744 1863959925 2677809493 33 d-19    /
      data bt02cs( 21) / -.1051301756 3758802637 9407414613 33 d-19    /
      data bt02cs( 22) / +.2182038199 1194813847 3010845013 33 d-20    /
      data bt02cs( 23) / -.4600470121 0362160577 2259054933 33 d-21    /
      data bt02cs( 24) / +.9840700692 5466818520 9536511999 99 d-22    /
      data bt02cs( 25) / -.2133403803 5728375844 7359863466 66 d-22    /
      data bt02cs( 26) / +.4683103642 3973365296 0662869333 33 d-23    /
      data bt02cs( 27) / -.1040021369 1985747236 5133823999 99 d-23    /
      data bt02cs( 28) / +.2334910567 7301510051 7777408000 00 d-24    /
      data bt02cs( 29) / -.5295682532 3318615788 0497493333 33 d-25    /
      data bt02cs( 30) / +.1212634195 2959756829 1962879999 99 d-25    /
      data bt02cs( 31) / -.2801889708 2289428760 2756266666 66 d-26    /
      data bt02cs( 32) / +.6529267898 7012873342 5937066666 66 d-27    /
      data bt02cs( 33) / -.1533798006 1873346427 8357333333 33 d-27    /
      data bt02cs( 34) / +.3630588430 6364536682 3594666666 66 d-28    /
      data bt02cs( 35) / -.8656075571 3629122479 1722666666 66 d-29    /
      data bt02cs( 36) / +.2077990997 2536284571 2383999999 99 d-29    /
      data bt02cs( 37) / -.5021117022 1417221674 3253333333 33 d-30    /
      data bt02cs( 38) / +.1220836027 9441714184 1919999999 99 d-30    /
      data bt02cs( 39) / -.2986005626 7039913454 2506666666 66 d-31    /
c
      data pi4 / 0.7853981633 9744830961 5660845819 876 d0 /
      data nbm0, nbt02, nbm02, nbth0, xmax / 4*0, 0.d0 /
c
      if (nbm0.ne.0) go to 10
      eta = 0.1*sngl(d1mach(3))
      nbm0 = initds (bm0cs, 37, eta)
      nbt02 = initds (bt02cs, 39, eta)
      nbm02 = initds (bm02cs, 40, eta)
      nbth0 = initds (bth0cs, 44, eta)
c
      xmax = 1.0d0/d1mach(4)
c
 10   if (x.lt.4.d0) call seteru (22hd9b0mp  x must be ge 4, 22, 1, 2)
c
      if (x.gt.8.d0) go to 20
      z = (128.d0/(x*x) - 5.d0)/3.d0
      ampl = (.75d0 + dcsevl (z, bm0cs, nbm0))/dsqrt(x)
      theta = x - pi4 + dcsevl (z, bt02cs, nbt02)/x
      return
c
 20   if (x.gt.xmax) call seteru (
     1  37hd9b0mp  no precision because x is big, 37, 2, 2)
c
      z = 128.d0/(x*x) - 1.d0
      ampl = (.75d0 + dcsevl (z, bm02cs, nbm02))/dsqrt(x)
      theta = x - pi4 + dcsevl (z, bth0cs, nbth0)/x
      return
c
      end