1 subroutine SBI0K0 (x, BI0, BK0, want, status) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5C>> 1998-10-29 SBI0K0 Krogh Moved external statement up for mangle. 6c>> 1996-04-27 SBI0K0 Krogh Changes to use .C. and C%%. 7C>> 1995-11-28 SBI0K0 Krogh Changes to simplify conversion to C. 8C>> 1994-10-20 SBI0K0 Krogh Changes to use M77CON 9C>> 1990-10-18 SBI0K0 WV Snyder JPL Use Cody K0 for small X 10C>> 1990-09-25 SBI0K0 WV Snyder JPL Use Fullerton codes from CMLIB 11C--S replaces "?": ?BI0K0, ?CSEVL, ?ERM1, ?INITS 12c 13c Compute hyperbolic Bessel functions I0 and K0. 14c Approximations for K0(X) on 0 < X < BOUNDK originally produced 15c by W. James Cody, ANL. Other approximations originally produced 16c by Wayne Fullerton, LASL. 17c 18c ***** Formal Arguments *********************************** 19c 20c X [in] is the argument at which the functions are to be evaluated. 21c BI0, BK0 [out] are the function values. 22c WANT [integer,in] indicates the functions to be computed, and their 23c scaling: 24c ABS(WANT) = 25c 1 means compute I0(X) 26c 2 means compute K0(X) 27c 3 means compute both of I0(X) and K0(X) 28c WANT < 0 means compute EXP(-X)*I0(X) and/or EXP(X)*K0(X). 29c WANT=0 or ABS(WANT) > 3 causes an error message to be produced. 30c STATUS [integer,out] indicates the outcome: 31c 0 means normal computation, 32c 1 means K0(X) is zero due to underflow, 33c < 0 means an error occurred: 34c -1 means WANT=0 or ABS(WANT)>3, 35c -2 means X was so big that I0(X) overflowed, 36c -3 means X was zero or negative and K0(X) is to be computed. 37c Negative values of STATUS are produced when the error message 38c processor is called with LEVEL=0; positive values of STATUS are 39c accompanied by LEVEL=-2. See the description of the error message 40c handler for a description of the error level effects. 41c If status = -2 then BI0 = the largest representable number; 42c if status = -3 then BK0 = the largest representable number. 43c ------------------------------------------------------------------ 44 real X, BI0, BK0 45 integer WANT, STATUS 46c 47c ***** External References ******************************** 48c 49 external SCSEVL, IERM1, SINITS, R1MACH, SERM1 50 real R1MACH, SCSEVL 51c 52c ***** Local Variables ************************************ 53c 54 real BI0CS(18), AI0CS(46), AI02CS(69) 55 real AK0CS(38), AK02CS(33) 56 real P(6), Q(2), PP(10), QQ(10), F(4), G(3) 57 real SUMF, SUMG, SUMP, SUMQ, Y, Z 58 real EXP10, EXPM10, LSQ2PI, LSQPI2 59 parameter (EXP10 = (+.2202646579480671651695790E+5)) 60C EXP10 = EXP(10) 61 parameter (EXPM10 = (+.4539992976248485153559152E-4)) 62C EXPM10 = EXP(-10) 63 parameter (LSQ2PI = (+.9189385332046727417803296E+0)) 64C LSQ2PI = LOG(SQRT (2 PI)) 65 parameter (LSQPI2 = (+.2257913526447274323630975E+0)) 66C LSQPI2 = LOG(SQRT(PI / 2)) 67 real BOUNDK, XIMAX, XIN, XISML, XKMAX, XKN, XKSML 68 integer I, NTI0, NTAI0, NTAI02, NTAK0, NTAK02 69 save NTI0, NTAI0, NTAI02, XIMAX, XISML 70 save BOUNDK, NTAK0, NTAK02, XKMAX, XKSML 71c 72c ***** Data *********************************************** 73c 74c SERIES FOR BI0 ON THE INTERVAL 0. TO 9.00000E+00 75c WITH WEIGHTED ERROR 9.51E-34 76c LOG WEIGHTED ERROR 33.02 77c SIGNIFICANT FIGURES REQUIRED 33.31 78c DECIMAL PLACES REQUIRED 33.65 79c 80 data BI0CS / -.7660547252839144951081894976243285E-1, 81 * +.1927337953993808269952408750881196E+01, 82 * +.2282644586920301338937029292330415E+00, 83 * +.1304891466707290428079334210691888E-01, 84 * +.4344270900816487451378682681026107E-03, 85 * +.9422657686001934663923171744118766E-05, 86 * +.1434006289510691079962091878179957E-06, 87 * +.1613849069661749069915419719994611E-08, 88 * +.1396650044535669699495092708142522E-10, 89 * +.9579451725505445344627523171893333E-13, 90 * +.5333981859862502131015107744000000E-15, 91 * +.2458716088437470774696785919999999E-17, 92 * +.9535680890248770026944341333333333E-20, 93 * +.3154382039721427336789333333333333E-22, 94 * +.9004564101094637431466666666666666E-25, 95 * +.2240647369123670016000000000000000E-27, 96 * +.4903034603242837333333333333333333E-30, 97 * +.9508172606122666666666666666666666E-33/ 98c 99c SERIES FOR AI0 ON THE INTERVAL 1.25000E-01 TO 3.33333E-01 100c WITH WEIGHTED ERROR 2.74E-32 101c LOG WEIGHTED ERROR 31.56 102c SIGNIFICANT FIGURES REQUIRED 30.15 103c DECIMAL PLACES REQUIRED 32.39 104c 105c++ Save data by elements if ~.C. 106 data AI0CS(1) / +.7575994494023795942729872037438E-1 / 107 data AI0CS(2) / +.7591380810823345507292978733204E-2 / 108 data AI0CS(3) / +.4153131338923750501863197491382E-3 / 109 data AI0CS(4) / +.1070076463439073073582429702170E-4 / 110 data AI0CS(5) / -.7901179979212894660750319485730E-5 / 111 data AI0CS(6) / -.7826143501438752269788989806909E-6 / 112 data AI0CS(7) / +.2783849942948870806381185389857E-6 / 113 data AI0CS(8) / +.8252472600612027191966829133198E-8 / 114 data AI0CS(9) / -.1204463945520199179054960891103E-7 / 115 data AI0CS(10) / +.1559648598506076443612287527928E-8 / 116 data AI0CS(11) / +.2292556367103316543477254802857E-9 / 117 data AI0CS(12) / -.1191622884279064603677774234478E-9 / 118 data AI0CS(13) / +.1757854916032409830218331247743E-10 / 119 data AI0CS(14) / +.1128224463218900517144411356824E-11 / 120 data AI0CS(15) / -.1146848625927298877729633876982E-11 / 121 data AI0CS(16) / +.2715592054803662872643651921606E-12 / 122 data AI0CS(17) / -.2415874666562687838442475720281E-13 / 123 data AI0CS(18) / -.6084469888255125064606099639224E-14 / 124 data AI0CS(19) / +.3145705077175477293708360267303E-14 / 125 data AI0CS(20) / -.7172212924871187717962175059176E-15 / 126 data AI0CS(21) / +.7874493403454103396083909603327E-16 / 127 data AI0CS(22) / +.1004802753009462402345244571839E-16 / 128 data AI0CS(23) / -.7566895365350534853428435888810E-17 / 129 data AI0CS(24) / +.2150380106876119887812051287845E-17 / 130 data AI0CS(25) / -.3754858341830874429151584452608E-18 / 131 data AI0CS(26) / +.2354065842226992576900757105322E-19 / 132 data AI0CS(27) / +.1114667612047928530226373355110E-19 / 133 data AI0CS(28) / -.5398891884396990378696779322709E-20 / 134 data AI0CS(29) / +.1439598792240752677042858404522E-20 / 135 data AI0CS(30) / -.2591916360111093406460818401962E-21 / 136 data AI0CS(31) / +.2238133183998583907434092298240E-22 / 137 data AI0CS(32) / +.5250672575364771172772216831999E-23 / 138 data AI0CS(33) / -.3249904138533230784173432285866E-23 / 139 data AI0CS(34) / +.9924214103205037927857284710400E-24 / 140 data AI0CS(35) / -.2164992254244669523146554299733E-24 / 141 data AI0CS(36) / +.3233609471943594083973332991999E-25 / 142 data AI0CS(37) / -.1184620207396742489824733866666E-26 / 143 data AI0CS(38) / -.1281671853950498650548338687999E-26 / 144 data AI0CS(39) / +.5827015182279390511605568853333E-27 / 145 data AI0CS(40) / -.1668222326026109719364501503999E-27 / 146 data AI0CS(41) / +.3625309510541569975700684800000E-28 / 147 data AI0CS(42) / -.5733627999055713589945958399999E-29 / 148 data AI0CS(43) / +.3736796722063098229642581333333E-30 / 149 data AI0CS(44) / +.1602073983156851963365512533333E-30 / 150 data AI0CS(45) / -.8700424864057229884522495999999E-31 / 151 data AI0CS(46) / +.2741320937937481145603413333333E-31 / 152c 153c SERIES FOR AI02 ON THE INTERVAL 0. TO 1.25000E-01 154c WITH WEIGHTED ERROR 1.97E-32 155c LOG WEIGHTED ERROR 31.71 156c SIGNIFICANT FIGURES REQUIRED 30.15 157c DECIMAL PLACES REQUIRED 32.63 158c 159c++ Save data by elements if ~.C. 160 data AI02CS(1) / +.5449041101410883160789609622680E-1 / 161 data AI02CS(2) / +.3369116478255694089897856629799E-2 / 162 data AI02CS(3) / +.6889758346916823984262639143011E-4 / 163 data AI02CS(4) / +.2891370520834756482966924023232E-5 / 164 data AI02CS(5) / +.2048918589469063741827605340931E-6 / 165 data AI02CS(6) / +.2266668990498178064593277431361E-7 / 166 data AI02CS(7) / +.3396232025708386345150843969523E-8 / 167 data AI02CS(8) / +.4940602388224969589104824497835E-9 / 168 data AI02CS(9) / +.1188914710784643834240845251963E-10 / 169 data AI02CS(10) / -.3149916527963241364538648629619E-10 / 170 data AI02CS(11) / -.1321581184044771311875407399267E-10 / 171 data AI02CS(12) / -.1794178531506806117779435740269E-11 / 172 data AI02CS(13) / +.7180124451383666233671064293469E-12 / 173 data AI02CS(14) / +.3852778382742142701140898017776E-12 / 174 data AI02CS(15) / +.1540086217521409826913258233397E-13 / 175 data AI02CS(16) / -.4150569347287222086626899720156E-13 / 176 data AI02CS(17) / -.9554846698828307648702144943125E-14 / 177 data AI02CS(18) / +.3811680669352622420746055355118E-14 / 178 data AI02CS(19) / +.1772560133056526383604932666758E-14 / 179 data AI02CS(20) / -.3425485619677219134619247903282E-15 / 180 data AI02CS(21) / -.2827623980516583484942055937594E-15 / 181 data AI02CS(22) / +.3461222867697461093097062508134E-16 / 182 data AI02CS(23) / +.4465621420296759999010420542843E-16 / 183 data AI02CS(24) / -.4830504485944182071255254037954E-17 / 184 data AI02CS(25) / -.7233180487874753954562272409245E-17 / 185 data AI02CS(26) / +.9921475412173698598880460939810E-18 / 186 data AI02CS(27) / +.1193650890845982085504399499242E-17 / 187 data AI02CS(28) / -.2488709837150807235720544916602E-18 / 188 data AI02CS(29) / -.1938426454160905928984697811326E-18 / 189 data AI02CS(30) / +.6444656697373443868783019493949E-19 / 190 data AI02CS(31) / +.2886051596289224326481713830734E-19 / 191 data AI02CS(32) / -.1601954907174971807061671562007E-19 / 192 data AI02CS(33) / -.3270815010592314720891935674859E-20 / 193 data AI02CS(34) / +.3686932283826409181146007239393E-20 / 194 data AI02CS(35) / +.1268297648030950153013595297109E-22 / 195 data AI02CS(36) / -.7549825019377273907696366644101E-21 / 196 data AI02CS(37) / +.1502133571377835349637127890534E-21 / 197 data AI02CS(38) / +.1265195883509648534932087992483E-21 / 198 data AI02CS(39) / -.6100998370083680708629408916002E-22 / 199 data AI02CS(40) / -.1268809629260128264368720959242E-22 / 200 data AI02CS(41) / +.1661016099890741457840384874905E-22 / 201 data AI02CS(42) / -.1585194335765885579379705048814E-23 / 202 data AI02CS(43) / -.3302645405968217800953817667556E-23 / 203 data AI02CS(44) / +.1313580902839239781740396231174E-23 / 204 data AI02CS(45) / +.3689040246671156793314256372804E-24 / 205 data AI02CS(46) / -.4210141910461689149219782472499E-24 / 206 data AI02CS(47) / +.4791954591082865780631714013730E-25 / 207 data AI02CS(48) / +.8459470390221821795299717074124E-25 / 208 data AI02CS(49) / -.4039800940872832493146079371810E-25 / 209 data AI02CS(50) / -.6434714653650431347301008504695E-26 / 210 data AI02CS(51) / +.1225743398875665990344647369905E-25 / 211 data AI02CS(52) / -.2934391316025708923198798211754E-26 / 212 data AI02CS(53) / -.1961311309194982926203712057289E-26 / 213 data AI02CS(54) / +.1503520374822193424162299003098E-26 / 214 data AI02CS(55) / -.9588720515744826552033863882069E-28 / 215 data AI02CS(56) / -.3483339380817045486394411085114E-27 / 216 data AI02CS(57) / +.1690903610263043673062449607256E-27 / 217 data AI02CS(58) / +.1982866538735603043894001157188E-28 / 218 data AI02CS(59) / -.5317498081491816214575830025284E-28 / 219 data AI02CS(60) / +.1803306629888392946235014503901E-28 / 220 data AI02CS(61) / +.6213093341454893175884053112422E-29 / 221 data AI02CS(62) / -.7692189292772161863200728066730E-29 / 222 data AI02CS(63) / +.1858252826111702542625560165963E-29 / 223 data AI02CS(64) / +.1237585142281395724899271545541E-29 / 224 data AI02CS(65) / -.1102259120409223803217794787792E-29 / 225 data AI02CS(66) / +.1886287118039704490077874479431E-30 / 226 data AI02CS(67) / +.2160196872243658913149031414060E-30 / 227 data AI02CS(68) / -.1605454124919743200584465949655E-30 / 228 data AI02CS(69) / +.1965352984594290603938848073318E-31 / 229c 230c ------------------------------------------------------------------- 231C 232C Coefficients for K0 from Cody for XSMALL .LE. ARG .LE. 1.0 233C 234c ------------------------------------------------------------------- 235 DATA P/ 5.8599221412826100000E-04, 1.3166052564989571850E-01, 236 1 1.1999463724910714109E+01, 4.6850901201934832188E+02, 237 2 5.9169059852270512312E+03, 2.4708152720399552679E+03/ 238 DATA Q/-2.4994418972832303646E+02, 2.1312714303849120380E+04/ 239 DATA F/-1.6414452837299064100E+00,-2.9601657892958843866E+02, 240 1 -1.7733784684952985886E+04,-4.0320340761145482298E+05/ 241 DATA G/-2.5064972445877992730E+02, 2.9865713163054025489E+04, 242 1 -1.6128136304458193998E+06/ 243c ------------------------------------------------------------------- 244C 245C Coefficients for K0 from Cody for 1.0 .LT. ARG .LT. BOUNDK 246C 247c ------------------------------------------------------------------- 248 DATA PP/ 1.1394980557384778174E+02, 3.6832589957340267940E+03, 249 1 3.1075408980684392399E+04, 1.0577068948034021957E+05, 250 2 1.7398867902565686251E+05, 1.5097646353289914539E+05, 251 3 7.1557062783764037541E+04, 1.8321525870183537725E+04, 252 4 2.3444738764199315021E+03, 1.1600249425076035558E+02/ 253 DATA QQ/ 2.0013443064949242491E+02, 4.4329628889746408858E+03, 254 1 3.1474655750295278825E+04, 9.7418829762268075784E+04, 255 2 1.5144644673520157801E+05, 1.2689839587977598727E+05, 256 3 5.8824616785857027752E+04, 1.4847228371802360957E+04, 257 4 1.8821890840982713696E+03, 9.2556599177304839811E+01/ 258c ------------------------------------------------------------------- 259C 260C SERIES FOR AK0 ON THE INTERVAL 1.25000E-01 TO 5.00000E-01 261C WITH WEIGHTED ERROR 2.85E-32 262C LOG WEIGHTED ERROR 31.54 263C SIGNIFICANT FIGURES REQUIRED 30.19 264C DECIMAL PLACES REQUIRED 32.33 265C 266c++ Save data by elements if ~.C. 267 data AK0CS(1) / -.7643947903327941424082978270088E-1 / 268 data AK0CS(2) / -.2235652605699819052023095550791E-1 / 269 data AK0CS(3) / +.7734181154693858235300618174047E-3 / 270 data AK0CS(4) / -.4281006688886099464452146435416E-4 / 271 data AK0CS(5) / +.3081700173862974743650014826660E-5 / 272 data AK0CS(6) / -.2639367222009664974067448892723E-6 / 273 data AK0CS(7) / +.2563713036403469206294088265742E-7 / 274 data AK0CS(8) / -.2742705549900201263857211915244E-8 / 275 data AK0CS(9) / +.3169429658097499592080832873403E-9 / 276 data AK0CS(10) / -.3902353286962184141601065717962E-10 / 277 data AK0CS(11) / +.5068040698188575402050092127286E-11 / 278 data AK0CS(12) / -.6889574741007870679541713557984E-12 / 279 data AK0CS(13) / +.9744978497825917691388201336831E-13 / 280 data AK0CS(14) / -.1427332841884548505389855340122E-13 / 281 data AK0CS(15) / +.2156412571021463039558062976527E-14 / 282 data AK0CS(16) / -.3349654255149562772188782058530E-15 / 283 data AK0CS(17) / +.5335260216952911692145280392601E-16 / 284 data AK0CS(18) / -.8693669980890753807639622378837E-17 / 285 data AK0CS(19) / +.1446404347862212227887763442346E-17 / 286 data AK0CS(20) / -.2452889825500129682404678751573E-18 / 287 data AK0CS(21) / +.4233754526232171572821706342400E-19 / 288 data AK0CS(22) / -.7427946526454464195695341294933E-20 / 289 data AK0CS(23) / +.1323150529392666866277967462400E-20 / 290 data AK0CS(24) / -.2390587164739649451335981465599E-21 / 291 data AK0CS(25) / +.4376827585923226140165712554666E-22 / 292 data AK0CS(26) / -.8113700607345118059339011413333E-23 / 293 data AK0CS(27) / +.1521819913832172958310378154666E-23 / 294 data AK0CS(28) / -.2886041941483397770235958613333E-24 / 295 data AK0CS(29) / +.5530620667054717979992610133333E-25 / 296 data AK0CS(30) / -.1070377329249898728591633066666E-25 / 297 data AK0CS(31) / +.2091086893142384300296328533333E-26 / 298 data AK0CS(32) / -.4121713723646203827410261333333E-27 / 299 data AK0CS(33) / +.8193483971121307640135680000000E-28 / 300 data AK0CS(34) / -.1642000275459297726780757333333E-28 / 301 data AK0CS(35) / +.3316143281480227195890346666666E-29 / 302 data AK0CS(36) / -.6746863644145295941085866666666E-30 / 303 data AK0CS(37) / +.1382429146318424677635413333333E-30 / 304 data AK0CS(38) / -.2851874167359832570811733333333E-31 / 305C 306C SERIES FOR AK02 ON THE INTERVAL 0. TO 1.25000E-01 307C WITH WEIGHTED ERROR 2.30E-32 308C LOG WEIGHTED ERROR 31.64 309C SIGNIFICANT FIGURES REQUIRED 29.68 310C DECIMAL PLACES REQUIRED 32.40 311C 312c++ Save data by elements if ~.C. 313 data AK02CS(1) / -.1201869826307592239839346212452E-1 / 314 data AK02CS(2) / -.9174852691025695310652561075713E-2 / 315 data AK02CS(3) / +.1444550931775005821048843878057E-3 / 316 data AK02CS(4) / -.4013614175435709728671021077879E-5 / 317 data AK02CS(5) / +.1567831810852310672590348990333E-6 / 318 data AK02CS(6) / -.7770110438521737710315799754460E-8 / 319 data AK02CS(7) / +.4611182576179717882533130529586E-9 / 320 data AK02CS(8) / -.3158592997860565770526665803309E-10 / 321 data AK02CS(9) / +.2435018039365041127835887814329E-11 / 322 data AK02CS(10) / -.2074331387398347897709853373506E-12 / 323 data AK02CS(11) / +.1925787280589917084742736504693E-13 / 324 data AK02CS(12) / -.1927554805838956103600347182218E-14 / 325 data AK02CS(13) / +.2062198029197818278285237869644E-15 / 326 data AK02CS(14) / -.2341685117579242402603640195071E-16 / 327 data AK02CS(15) / +.2805902810643042246815178828458E-17 / 328 data AK02CS(16) / -.3530507631161807945815482463573E-18 / 329 data AK02CS(17) / +.4645295422935108267424216337066E-19 / 330 data AK02CS(18) / -.6368625941344266473922053461333E-20 / 331 data AK02CS(19) / +.9069521310986515567622348800000E-21 / 332 data AK02CS(20) / -.1337974785423690739845005311999E-21 / 333 data AK02CS(21) / +.2039836021859952315522088960000E-22 / 334 data AK02CS(22) / -.3207027481367840500060869973333E-23 / 335 data AK02CS(23) / +.5189744413662309963626359466666E-24 / 336 data AK02CS(24) / -.8629501497540572192964607999999E-25 / 337 data AK02CS(25) / +.1472161183102559855208038400000E-25 / 338 data AK02CS(26) / -.2573069023867011283812351999999E-26 / 339 data AK02CS(27) / +.4601774086643516587376640000000E-27 / 340 data AK02CS(28) / -.8411555324201093737130666666666E-28 / 341 data AK02CS(29) / +.1569806306635368939301546666666E-28 / 342 data AK02CS(30) / -.2988226453005757788979199999999E-29 / 343 data AK02CS(31) / +.5796831375216836520618666666666E-30 / 344 data AK02CS(32) / -.1145035994347681332155733333333E-30 / 345 data AK02CS(33) / +.2301266594249682802005333333333E-31 / 346c 347 data NTI0 /0/ 348c 349c ***** Statement Functions ******************************** 350c 351 xin(z) = ximax + 0.5*log(z/(1.0+1.0/(8.0*z))**2) 352 xkn(z) = xkmax - 0.5*log(z/(1.0-1.0/(8.0*z))**2) 353c 354c ***** Executable Statements ****************************** 355c 356 if (nti0 .le. 0) then 357 z = 0.1*R1MACH(3) 358 boundk = 1.757 - 6.22e-3 * log(z) 359 call SINITS (bi0cs, 18, z, nti0 ) 360 call SINITS (ai0cs, 46, z, ntai0 ) 361 call SINITS (ai02cs, 69, z, ntai02) 362 call SINITS (ak0cs, 38, z, ntak0 ) 363 call SINITS (ak02cs, 33, z, ntak02) 364 xisml = sqrt (80.0E0*z) 365 ximax = log(R1MACH(2)) + lsq2pi 366 ximax = xin(xin(ximax)) 367 xksml = sqrt (10.0*z*p(6)/p(5)) 368c xksml = sqrt (10.0*z*min(p(6)/p(5),q(2)/q(1))) 369 xkmax = lsqpi2 - log(R1MACH(1)) 370 xkmax = xkn(xkn(xkn(xkmax))) 371 end if 372c 373 if (want.eq.0 .or. abs(want).gt.3) then 374 call ierm1('SBI0K0',-1,0,'WANT = 0 OR ABS(WANT) > 3','WANT', 375 1 want,'.') 376 status=-1 377 return 378 end if 379 status = 0 380c 381c Compute I0 if requested. 382c 383 if (abs(want) .ne. 2) then 384 y = abs(x) 385 if (y.le.3.0) then 386 if (y.le.xisml) then 387 BI0 = 1.0E0 388 else 389 BI0 = 2.75E0 + SCSEVL(y*y/4.5E0-1.0E0, bi0cs, nti0) 390 end if 391 if (want.lt.0) BI0 = exp(-y) * BI0 392 else 393 if (y .le. 8.0E0) then 394 BI0 = SCSEVL ((48.0E0/y-11.0E0)/5.0E0, ai0cs, ntai0) 395 else 396 BI0 = SCSEVL (16.0E0/y-1.0E0, ai02cs, ntai02) 397 end if 398 BI0 = (0.375E0 + BI0) / sqrt(y) 399 if (want .gt. 0) then 400 if (y .gt. ximax) then 401 call SERM1('SBI0K0',-2,0,'ABS(X) SO BIG I0 OVERFLOWS', 402 1 'X',x,'.') 403 status = -2 404 BI0 = R1MACH(2) 405c y > ximax => y > xkmax 406 BK0 = 0.0 407 if (x .gt. 0.0) return 408 else 409 BI0 = (exp(y-10.0) * BI0) * exp10 410 end if 411 end if 412 end if 413 end if 414c 415c Compute K0 if requested. 416c 417 if (abs(want) .gt. 1) then 418 if (x .le. 0.0E0) then 419 call SERM1('SBI0K0',-3,0,'X IS ZERO OR NEGATIVE','X',x,'.') 420 status = -3 421 BK0 = R1MACH(2) 422 else if (x .le. 1.0) then 423c 0.0 < x <= 1.0 424 BK0 = log(x) 425 if (x .lt. xksml) then 426c 0.0 <= x < xksml 427 BK0 = p(6)/q(2) - BK0 428 else 429c xksml <= x <= 1.0 430 y = x*x 431 sump = ((((p(1)*y+p(2))*y+p(3))*y+p(4))*y+p(5))*y+p(6) 432 sumq = (y + q(1))*y + q(2) 433 sumf = ((f(1)*y + f(2))*y + f(3))*y + f(4) 434 sumg = ((y + g(1))*y + g(2))*y + g(3) 435 BK0 = sump/sumq - y * sumf * BK0 / sumg - BK0 436 end if 437 if (want .lt. 0) BK0 = exp(x) * BK0 438 else if (x .le. boundk) then 439c 1.0 < x <= boundk 440 y = 1.0 / x 441 sump = pp(1) 442 do 120 i = 2, 10 443 sump = sump*y + pp(i) 444120 continue 445 sumq = y 446 do 140 i = 1, 9 447 sumq = (sumq + qq(i))*y 448140 continue 449 sumq = sumq + qq(10) 450 BK0 = sump / (sumq * sqrt(x)) 451 if (want .gt. 0) BK0 = exp(-x) * BK0 452 else 453c boundk < x 454 y = 16.0E0 / x 455 if (x .le. 8.0E0) then 456 BK0 = SCSEVL ((y-5.0E0)/3.0E0, ak0cs, ntak0) 457 else 458 BK0 = SCSEVL (y-1.0E0, ak02cs, ntak02) 459 end if 460 BK0 = (1.25E0 + BK0) / sqrt(x) 461 if (want .gt. 0) then 462 if (x .gt. xkmax) then 463 call SERM1('SBI0K0',1,-2,'X SO BIG K0 UNDERFLOWS','X', 464 1 x,'.') 465 BK0 = 0.0 466 status = 1 467 else 468 BK0 = (exp(10.0-x) * BK0) * expm10 469 end if 470 end if 471 end if 472 end if 473C 474 return 475c 476 end 477