subroutine besscheb(betal,ki,kii) real*8 betal, ki, kii c INPUT: betal = Log10(kT/m c**2) c Returns ki = exp(z)*KI, kii = exp(z)*KII, c calculated with an expansion in Chebyshev polynomials. c KI and KII are modified Bessel functions of the second kind. c We use the expansion given by Tooper (1969, ApJ 156, 1075) c for z>8 and the expansion given by Clenshaw (1962, Mathematical c Tables, vol. 5, Chebyshev Series for Mathematical Functions, c London: Her Majesty's Stationery Office. 510 N279m v.5) for z<8. c Note: There is a mistake in Tooper's (13)-(14) Tables c The formula for K should be w/o the factor pi. c Written by Juan Antonio Miralles, July, 1995 real*8 ak0(0:14), aki(0:14) real*8 bi0(0:36), bk0(0:36), bii(0:34), bki(0:34) real*8 k0, i0, ii real*8 beta, z, x, y, bkp1, bkp2, bkp3, bk integer k logical initialized data initialized /.false./ save ak0, aki ,bi0, bk0, bii, bki save initialized if (.not.initialized) then ak0(0) = 2.487981301736924078 ak0(1) = -9.17485269101569531d-3 ak0(2) = 1.4445509317750058d-4 ak0(3) = -4.01361417543571d-6 ak0(4) = 1.5678318108523d-7 ak0(5) = -7.77011043852d-9 ak0(6) = 4.6111825762d-10 ak0(7) = -3.158592998d-11 ak0(8) = 2.43501804d-12 ak0(9) = -2.0743314d-13 ak0(10) = 1.925787d-14 ak0(11) = -1.92755d-15 ak0(12) = 2.0622d-16 ak0(13) = -2.342d-17 ak0(14) = 2.81d-18 aki(0) = 2.563793083437390010 aki(1) = 2.832887813049720936d-2 aki(2) = -2.4753706739052503d-4 aki(3) = 5.77197245160725d-6 aki(4) = -2.0689392195365d-7 aki(5) = 9.73998344138d-9 aki(6) = -5.5853361404d-10 aki(7) = 3.732996634d-11 aki(8) = -2.82505196d-12 aki(9) = 2.3720190d-13 aki(10) = -2.176677d-14 aki(11) = 2.15791d-15 aki(12) = -2.2902d-16 aki(13) = 2.583d-17 aki(14) = -3.08d-18 bi0(0) = 255.46687962436216712602 bi0(2) = 190.49432017274284419322 bi0(4) = 82.48903274402409961321 bi0(6) = 22.27481924246223087742 bi0(8) = 4.01167376017934853351 bi0(10) = 0.50949336543998287079 bi0(12) = 0.04771874879817413524 bi0(14) = 0.00341633176601234095 bi0(16) = 0.00019246935968811366 bi0(18) = 0.00000873831549662236 bi0(20) = 0.00000032609105057896 bi0(22) = 0.00000001016972672769 bi0(24) = 0.00000000026882812895 bi0(26) = 0.00000000000609689280 bi0(28) = 0.00000000000011989083 bi0(30) = 0.00000000000000206305 bi0(32) = 0.00000000000000003132 bi0(34) = 0.00000000000000000042 bi0(36) = 0.00000000000000000001 bk0(0) = -21.05766017740244024847 bk0(2) = - 4.56343358644839500735 bk0(4) = 8.00536886870033477053 bk0(6) = 5.28363286687392000748 bk0(8) = 1.51153567602922791353 bk0(10) = 0.25908443243490019747 bk0(12) = 0.03008072242051187452 bk0(14) = 0.00253630818808619901 bk0(16) = 0.00016270837904302330 bk0(18) = 0.00000821602593993066 bk0(20) = 0.00000033519525563133 bk0(22) = 0.00000001128121138760 bk0(24) = 0.00000000031858797963 bk0(26) = 0.00000000000765757438 bk0(28) = 0.00000000000015855413 bk0(30) = 0.00000000000000285752 bk0(32) = 0.00000000000000004523 bk0(34) = 0.00000000000000000063 bk0(36) = 0.00000000000000000001 bii(0) = 259.89023780647729173120 bii(2) = 181.31261604057026539103 bii(4) = 69.39591763373444753798 bii(6) = 16.33455055252206616461 bii(8) = 2.57145990634775490572 bii(10) = 0.28785551180467203057 bii(12) = 0.02399307914784055176 bii(14) = 0.00154301901562721914 bii(16) = 0.00007875678575416515 bii(18) = 0.00000326413812230986 bii(20) = 0.00000011194628456389 bii(22) = 0.00000000322761652023 bii(24) = 0.00000000007929055929 bii(26) = 0.00000000000167897282 bii(28) = 0.00000000000003095296 bii(30) = 0.00000000000000050120 bii(32) = 0.00000000000000000718 bii(34) = 0.00000000000000000009 bki(0) = -26.68809548086266779437 bki(2) = - 1.83923922428619943219 bki(4) = 9.36161783139538867938 bki(6) = 4.66638702686284185037 bki(8) = 1.10146199300485220124 bki(10) = 0.16107430165614782390 bki(12) = 0.01630004928981641757 bki(14) = 0.00121705699451574089 bki(16) = 0.00007001062785475753 bki(18) = 0.00000320251069193505 bki(20) = 0.00000011936797074664 bki(22) = 0.00000000369678327836 bki(24) = 0.00000000009665975200 bki(26) = 0.00000000000216255319 bki(28) = 0.00000000000004187279 bki(30) = 0.00000000000000070860 bki(32) = 0.00000000000000001057 bki(34) = 0.00000000000000000014 initialized = .true. endif beta = 10.**betal z = 1./beta x = 16./z-1. y = z/8. if (z.gt.8.)then c Tooper's expansion bkp2 = 0. bkp1 = 0. do k = 14,0,-1 bkp3 = bkp2 bk = 2.*x*bkp1-bkp2+ak0(k) bkp2 = bkp1 bkp1 = bk enddo k0 = 0.5*(bkp1 - bkp3) bkp2 = 0. bkp1 = 0. do k = 14,0,-1 bkp3 = bkp2 bk = 2.*x*bkp1-bkp2+aki(k) bkp2 = bkp1 bkp1 = bk enddo ki = 0.5*(bkp1 - bkp3) kii = k0+2*ki/z c Note the following 2 do NOT contain factor exp(-z) ki = ki/sqrt(z) kii = kii/sqrt(z) else c Clenshaw's expansion bkp2 = 0. bkp1 = 0. do k = 18,0,-1 bkp3 = bkp2 bk = 2.*(2.*y*y-1.)*bkp1-bkp2+bi0(2*k) bkp2 = bkp1 bkp1 = bk enddo i0 = 0.5*(bkp1 - bkp3) bkp2 = 0. bkp1 = 0. do k = 18,0,-1 bkp3 = bkp2 bk = 2.*(2.*y*y-1.)*bkp1-bkp2+bk0(2*k) bkp2 = bkp1 bkp1 = bk enddo k0 = 0.5*(bkp1 - bkp3) k0 = -log(y)*i0 + k0 bkp2 = 0. bkp1 = 0. do k = 17,0,-1 bkp3 = bkp2 bk = 2.*(2.*y*y-1.)*bkp1-bkp2+bii(2*k) bkp2 = bkp1 bkp1 = bk enddo ii = 0.5*(bkp1 - bkp3) ii = z*ii/8. bkp2 = 0. bkp1 = 0. do k = 17,0,-1 bkp3 = bkp2 bk = 2.*(2.*y*y-1.)*bkp1-bkp2+bki(2*k) bkp2 = bkp1 bkp1 = bk enddo ki = 0.5*(bkp1 - bkp3) ki = log(y)*ii + 1./z - z*ki/8. kii = k0+2*ki/z c Note the following 2 do NOT contain factor exp(-z) ki = ki*exp(z) kii = kii*exp(z) endif return end