/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:57 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sbi0k0.h" #include #include /* PARAMETER translations */ #define EXP10 (.2202646579480671651695790e5) #define EXPM10 (.4539992976248485153559152e-4) #define LSQ2PI (.9189385332046727417803296e0) #define LSQPI2 (.2257913526447274323630975e0) /* end of PARAMETER translations */ void /*FUNCTION*/ sbi0k0( float x, float *bi0, float *bk0, long want, long *status) { long int _l0, i; static long int ntai0, ntai02, ntak0, ntak02; float sumf, sumg, sump, sumq, y, z; static float boundk, ximax, xisml, xkmax, xksml; static float bi0cs[18]={-.7660547252839144951081894976243285e-1, .1927337953993808269952408750881196e01,.2282644586920301338937029292330415e00, .1304891466707290428079334210691888e-01,.4344270900816487451378682681026107e-03, .9422657686001934663923171744118766e-05,.1434006289510691079962091878179957e-06, .1613849069661749069915419719994611e-08,.1396650044535669699495092708142522e-10, .9579451725505445344627523171893333e-13,.5333981859862502131015107744000000e-15, .2458716088437470774696785919999999e-17,.9535680890248770026944341333333333e-20, .3154382039721427336789333333333333e-22,.9004564101094637431466666666666666e-25, .2240647369123670016000000000000000e-27,.4903034603242837333333333333333333e-30, .9508172606122666666666666666666666e-33}; static float ai0cs[46]={.7575994494023795942729872037438e-1,.7591380810823345507292978733204e-2, .4153131338923750501863197491382e-3,.1070076463439073073582429702170e-4, -.7901179979212894660750319485730e-5,-.7826143501438752269788989806909e-6, .2783849942948870806381185389857e-6,.8252472600612027191966829133198e-8, -.1204463945520199179054960891103e-7,.1559648598506076443612287527928e-8, .2292556367103316543477254802857e-9,-.1191622884279064603677774234478e-9, .1757854916032409830218331247743e-10,.1128224463218900517144411356824e-11, -.1146848625927298877729633876982e-11,.2715592054803662872643651921606e-12, -.2415874666562687838442475720281e-13,-.6084469888255125064606099639224e-14, .3145705077175477293708360267303e-14,-.7172212924871187717962175059176e-15, .7874493403454103396083909603327e-16,.1004802753009462402345244571839e-16, -.7566895365350534853428435888810e-17,.2150380106876119887812051287845e-17, -.3754858341830874429151584452608e-18,.2354065842226992576900757105322e-19, .1114667612047928530226373355110e-19,-.5398891884396990378696779322709e-20, .1439598792240752677042858404522e-20,-.2591916360111093406460818401962e-21, .2238133183998583907434092298240e-22,.5250672575364771172772216831999e-23, -.3249904138533230784173432285866e-23,.9924214103205037927857284710400e-24, -.2164992254244669523146554299733e-24,.3233609471943594083973332991999e-25, -.1184620207396742489824733866666e-26,-.1281671853950498650548338687999e-26, .5827015182279390511605568853333e-27,-.1668222326026109719364501503999e-27, .3625309510541569975700684800000e-28,-.5733627999055713589945958399999e-29, .3736796722063098229642581333333e-30,.1602073983156851963365512533333e-30, -.8700424864057229884522495999999e-31,.2741320937937481145603413333333e-31}; static float ai02cs[69]={.5449041101410883160789609622680e-1,.3369116478255694089897856629799e-2, .6889758346916823984262639143011e-4,.2891370520834756482966924023232e-5, .2048918589469063741827605340931e-6,.2266668990498178064593277431361e-7, .3396232025708386345150843969523e-8,.4940602388224969589104824497835e-9, .1188914710784643834240845251963e-10,-.3149916527963241364538648629619e-10, -.1321581184044771311875407399267e-10,-.1794178531506806117779435740269e-11, .7180124451383666233671064293469e-12,.3852778382742142701140898017776e-12, .1540086217521409826913258233397e-13,-.4150569347287222086626899720156e-13, -.9554846698828307648702144943125e-14,.3811680669352622420746055355118e-14, .1772560133056526383604932666758e-14,-.3425485619677219134619247903282e-15, -.2827623980516583484942055937594e-15,.3461222867697461093097062508134e-16, .4465621420296759999010420542843e-16,-.4830504485944182071255254037954e-17, -.7233180487874753954562272409245e-17,.9921475412173698598880460939810e-18, .1193650890845982085504399499242e-17,-.2488709837150807235720544916602e-18, -.1938426454160905928984697811326e-18,.6444656697373443868783019493949e-19, .2886051596289224326481713830734e-19,-.1601954907174971807061671562007e-19, -.3270815010592314720891935674859e-20,.3686932283826409181146007239393e-20, .1268297648030950153013595297109e-22,-.7549825019377273907696366644101e-21, .1502133571377835349637127890534e-21,.1265195883509648534932087992483e-21, -.6100998370083680708629408916002e-22,-.1268809629260128264368720959242e-22, .1661016099890741457840384874905e-22,-.1585194335765885579379705048814e-23, -.3302645405968217800953817667556e-23,.1313580902839239781740396231174e-23, .3689040246671156793314256372804e-24,-.4210141910461689149219782472499e-24, .4791954591082865780631714013730e-25,.8459470390221821795299717074124e-25, -.4039800940872832493146079371810e-25,-.6434714653650431347301008504695e-26, .1225743398875665990344647369905e-25,-.2934391316025708923198798211754e-26, -.1961311309194982926203712057289e-26,.1503520374822193424162299003098e-26, -.9588720515744826552033863882069e-28,-.3483339380817045486394411085114e-27, .1690903610263043673062449607256e-27,.1982866538735603043894001157188e-28, -.5317498081491816214575830025284e-28,.1803306629888392946235014503901e-28, .6213093341454893175884053112422e-29,-.7692189292772161863200728066730e-29, .1858252826111702542625560165963e-29,.1237585142281395724899271545541e-29, -.1102259120409223803217794787792e-29,.1886287118039704490077874479431e-30, .2160196872243658913149031414060e-30,-.1605454124919743200584465949655e-30, .1965352984594290603938848073318e-31}; static float p[6]={5.8599221412826100000e-04,1.3166052564989571850e-01, 1.1999463724910714109e01,4.6850901201934832188e02,5.9169059852270512312e03, 2.4708152720399552679e03}; static float q[2]={-2.4994418972832303646e02,2.1312714303849120380e04}; static float f[4]={-1.6414452837299064100e00,-2.9601657892958843866e02, -1.7733784684952985886e04,-4.0320340761145482298e05}; static float g[3]={-2.5064972445877992730e02,2.9865713163054025489e04, -1.6128136304458193998e06}; static float pp[10]={1.1394980557384778174e02,3.6832589957340267940e03, 3.1075408980684392399e04,1.0577068948034021957e05,1.7398867902565686251e05, 1.5097646353289914539e05,7.1557062783764037541e04,1.8321525870183537725e04, 2.3444738764199315021e03,1.1600249425076035558e02}; static float qq[10]={2.0013443064949242491e02,4.4329628889746408858e03, 3.1474655750295278825e04,9.7418829762268075784e04,1.5144644673520157801e05, 1.2689839587977598727e05,5.8824616785857027752e04,1.4847228371802360957e04, 1.8821890840982713696e03,9.2556599177304839811e01}; static float ak0cs[38]={-.7643947903327941424082978270088e-1,-.2235652605699819052023095550791e-1, .7734181154693858235300618174047e-3,-.4281006688886099464452146435416e-4, .3081700173862974743650014826660e-5,-.2639367222009664974067448892723e-6, .2563713036403469206294088265742e-7,-.2742705549900201263857211915244e-8, .3169429658097499592080832873403e-9,-.3902353286962184141601065717962e-10, .5068040698188575402050092127286e-11,-.6889574741007870679541713557984e-12, .9744978497825917691388201336831e-13,-.1427332841884548505389855340122e-13, .2156412571021463039558062976527e-14,-.3349654255149562772188782058530e-15, .5335260216952911692145280392601e-16,-.8693669980890753807639622378837e-17, .1446404347862212227887763442346e-17,-.2452889825500129682404678751573e-18, .4233754526232171572821706342400e-19,-.7427946526454464195695341294933e-20, .1323150529392666866277967462400e-20,-.2390587164739649451335981465599e-21, .4376827585923226140165712554666e-22,-.8113700607345118059339011413333e-23, .1521819913832172958310378154666e-23,-.2886041941483397770235958613333e-24, .5530620667054717979992610133333e-25,-.1070377329249898728591633066666e-25, .2091086893142384300296328533333e-26,-.4121713723646203827410261333333e-27, .8193483971121307640135680000000e-28,-.1642000275459297726780757333333e-28, .3316143281480227195890346666666e-29,-.6746863644145295941085866666666e-30, .1382429146318424677635413333333e-30,-.2851874167359832570811733333333e-31}; static float ak02cs[33]={-.1201869826307592239839346212452e-1, -.9174852691025695310652561075713e-2,.1444550931775005821048843878057e-3, -.4013614175435709728671021077879e-5,.1567831810852310672590348990333e-6, -.7770110438521737710315799754460e-8,.4611182576179717882533130529586e-9, -.3158592997860565770526665803309e-10,.2435018039365041127835887814329e-11, -.2074331387398347897709853373506e-12,.1925787280589917084742736504693e-13, -.1927554805838956103600347182218e-14,.2062198029197818278285237869644e-15, -.2341685117579242402603640195071e-16,.2805902810643042246815178828458e-17, -.3530507631161807945815482463573e-18,.4645295422935108267424216337066e-19, -.6368625941344266473922053461333e-20,.9069521310986515567622348800000e-21, -.1337974785423690739845005311999e-21,.2039836021859952315522088960000e-22, -.3207027481367840500060869973333e-23,.5189744413662309963626359466666e-24, -.8629501497540572192964607999999e-25,.1472161183102559855208038400000e-25, -.2573069023867011283812351999999e-26,.4601774086643516587376640000000e-27, -.8411555324201093737130666666666e-28,.1569806306635368939301546666666e-28, -.2988226453005757788979199999999e-29,.5796831375216836520618666666666e-30, -.1145035994347681332155733333333e-30,.2301266594249682802005333333333e-31}; static long nti0 = 0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Ai02cs = &ai02cs[0] - 1; float *const Ai0cs = &ai0cs[0] - 1; float *const Ak02cs = &ak02cs[0] - 1; float *const Ak0cs = &ak0cs[0] - 1; float *const Bi0cs = &bi0cs[0] - 1; float *const F = &f[0] - 1; float *const G = &g[0] - 1; float *const P = &p[0] - 1; float *const Pp = &pp[0] - 1; float *const Q = &q[0] - 1; float *const Qq = &qq[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1998-10-29 SBI0K0 Krogh Moved external statement up for mangle. *>> 1996-04-27 SBI0K0 Krogh Changes to use .C. and C%%. *>> 1995-11-28 SBI0K0 Krogh Changes to simplify conversion to C. *>> 1994-10-20 SBI0K0 Krogh Changes to use M77CON *>> 1990-10-18 SBI0K0 WV Snyder JPL Use Cody K0 for small X *>> 1990-09-25 SBI0K0 WV Snyder JPL Use Fullerton codes from CMLIB *--S replaces "?": ?BI0K0, ?CSEVL, ?ERM1, ?INITS * * Compute hyperbolic Bessel functions I0 and K0. * Approximations for K0(X) on 0 < X < BOUNDK originally produced * by W. James Cody, ANL. Other approximations originally produced * by Wayne Fullerton, LASL. * * ***** Formal Arguments *********************************** * * X [in] is the argument at which the functions are to be evaluated. * BI0, BK0 [out] are the function values. * WANT [integer,in] indicates the functions to be computed, and their * scaling: * ABS(WANT) = * 1 means compute I0(X) * 2 means compute K0(X) * 3 means compute both of I0(X) and K0(X) * WANT < 0 means compute EXP(-X)*I0(X) and/or EXP(X)*K0(X). * WANT=0 or ABS(WANT) > 3 causes an error message to be produced. * STATUS [integer,out] indicates the outcome: * 0 means normal computation, * 1 means K0(X) is zero due to underflow, * < 0 means an error occurred: * -1 means WANT=0 or ABS(WANT)>3, * -2 means X was so big that I0(X) overflowed, * -3 means X was zero or negative and K0(X) is to be computed. * Negative values of STATUS are produced when the error message * processor is called with LEVEL=0; positive values of STATUS are * accompanied by LEVEL=-2. See the description of the error message * handler for a description of the error level effects. * If status = -2 then BI0 = the largest representable number; * if status = -3 then BK0 = the largest representable number. * ------------------------------------------------------------------ */ /* ***** External References ******************************** * */ /* ***** Local Variables ************************************ * */ /* EXP10 = EXP(10) */ /* EXPM10 = EXP(-10) */ /* LSQ2PI = LOG(SQRT (2 PI)) */ /* LSQPI2 = LOG(SQRT(PI / 2)) */ /* ***** Data *********************************************** * * SERIES FOR BI0 ON THE INTERVAL 0. TO 9.00000E+00 * WITH WEIGHTED ERROR 9.51E-34 * LOG WEIGHTED ERROR 33.02 * SIGNIFICANT FIGURES REQUIRED 33.31 * DECIMAL PLACES REQUIRED 33.65 * */ /* SERIES FOR AI0 ON THE INTERVAL 1.25000E-01 TO 3.33333E-01 * WITH WEIGHTED ERROR 2.74E-32 * LOG WEIGHTED ERROR 31.56 * SIGNIFICANT FIGURES REQUIRED 30.15 * DECIMAL PLACES REQUIRED 32.39 * *++ Save data by elements if ~.C. */ /* SERIES FOR AI02 ON THE INTERVAL 0. TO 1.25000E-01 * WITH WEIGHTED ERROR 1.97E-32 * LOG WEIGHTED ERROR 31.71 * SIGNIFICANT FIGURES REQUIRED 30.15 * DECIMAL PLACES REQUIRED 32.63 * *++ Save data by elements if ~.C. */ /* ------------------------------------------------------------------- * * Coefficients for K0 from Cody for XSMALL .LE. ARG .LE. 1.0 * * ------------------------------------------------------------------- */ /* ------------------------------------------------------------------- * * Coefficients for K0 from Cody for 1.0 .LT. ARG .LT. BOUNDK * * ------------------------------------------------------------------- */ /* ------------------------------------------------------------------- * * SERIES FOR AK0 ON THE INTERVAL 1.25000E-01 TO 5.00000E-01 * WITH WEIGHTED ERROR 2.85E-32 * LOG WEIGHTED ERROR 31.54 * SIGNIFICANT FIGURES REQUIRED 30.19 * DECIMAL PLACES REQUIRED 32.33 * *++ Save data by elements if ~.C. */ /* SERIES FOR AK02 ON THE INTERVAL 0. TO 1.25000E-01 * WITH WEIGHTED ERROR 2.30E-32 * LOG WEIGHTED ERROR 31.64 * SIGNIFICANT FIGURES REQUIRED 29.68 * DECIMAL PLACES REQUIRED 32.40 * *++ Save data by elements if ~.C. */ /* ***** Statement Functions ******************************** * */ #define XIN(z) ((float)(ximax + 0.5*logf( (z)/powif(1.0 + 1.0/(8.0*\ (z)),2) ))) #define XKN(z) ((float)(xkmax - 0.5*logf( (z)/powif(1.0 - 1.0/(8.0*\ (z)),2) ))) /* ***** Executable Statements ****************************** * */ if (nti0 <= 0) { z = 0.1*FLT_EPSILON/FLT_RADIX; boundk = 1.757 - 6.22e-3*logf( z ); sinits( bi0cs, 18, z, &nti0 ); sinits( ai0cs, 46, z, &ntai0 ); sinits( ai02cs, 69, z, &ntai02 ); sinits( ak0cs, 38, z, &ntak0 ); sinits( ak02cs, 33, z, &ntak02 ); xisml = sqrtf( 80.0e0*z ); ximax = logf( FLT_MAX ) + LSQ2PI; ximax = XIN( XIN( ximax ) ); xksml = sqrtf( 10.0*z*P[6]/P[5] ); /* xksml = sqrt (10.0*z*min(p(6)/p(5),q(2)/q(1))) */ xkmax = LSQPI2 - logf( FLT_MIN ); xkmax = XKN( XKN( XKN( xkmax ) ) ); } if (want == 0 || labs( want ) > 3) { ierm1( "SBI0K0", -1, 0, "WANT = 0 OR ABS(WANT) > 3", "WANT" , want, '.' ); *status = -1; return; } *status = 0; /* Compute I0 if requested. * */ if (labs( want ) != 2) { y = fabsf( x ); if (y <= 3.0) { if (y <= xisml) { *bi0 = 1.0e0; } else { *bi0 = 2.75e0 + scsevl( y*y/4.5e0 - 1.0e0, bi0cs, nti0 ); } if (want < 0) *bi0 *= expf( -y ); } else { if (y <= 8.0e0) { *bi0 = scsevl( (48.0e0/y - 11.0e0)/5.0e0, ai0cs, ntai0 ); } else { *bi0 = scsevl( 16.0e0/y - 1.0e0, ai02cs, ntai02 ); } *bi0 = (0.375e0 + *bi0)/sqrtf( y ); if (want > 0) { if (y > ximax) { serm1( "SBI0K0", -2, 0, "ABS(X) SO BIG I0 OVERFLOWS" , "X", x, '.' ); *status = -2; *bi0 = FLT_MAX; /* y > ximax => y > xkmax */ *bk0 = 0.0; if (x > 0.0) return; } else { *bi0 = (expf( y - 10.0 )**bi0)*EXP10; } } } } /* Compute K0 if requested. * */ if (labs( want ) > 1) { if (x <= 0.0e0) { serm1( "SBI0K0", -3, 0, "X IS ZERO OR NEGATIVE", "X", x, '.' ); *status = -3; *bk0 = FLT_MAX; } else if (x <= 1.0) { /* 0.0 < x <= 1.0 */ *bk0 = logf( x ); if (x < xksml) { /* 0.0 <= x < xksml */ *bk0 = P[6]/Q[2] - *bk0; } else { /* xksml <= x <= 1.0 */ y = x*x; sump = ((((P[1]*y + P[2])*y + P[3])*y + P[4])*y + P[5])*y + P[6]; sumq = (y + Q[1])*y + Q[2]; sumf = ((F[1]*y + F[2])*y + F[3])*y + F[4]; sumg = ((y + G[1])*y + G[2])*y + G[3]; *bk0 = sump/sumq - y*sumf**bk0/sumg - *bk0; } if (want < 0) *bk0 *= expf( x ); } else if (x <= boundk) { /* 1.0 < x <= boundk */ y = 1.0/x; sump = Pp[1]; for (i = 2; i <= 10; i++) { sump = sump*y + Pp[i]; } sumq = y; for (i = 1; i <= 9; i++) { sumq = (sumq + Qq[i])*y; } sumq += Qq[10]; *bk0 = sump/(sumq*sqrtf( x )); if (want > 0) *bk0 *= expf( -x ); } else { /* boundk < x */ y = 16.0e0/x; if (x <= 8.0e0) { *bk0 = scsevl( (y - 5.0e0)/3.0e0, ak0cs, ntak0 ); } else { *bk0 = scsevl( y - 1.0e0, ak02cs, ntak02 ); } *bk0 = (1.25e0 + *bk0)/sqrtf( x ); if (want > 0) { if (x > xkmax) { serm1( "SBI0K0", 1, -2, "X SO BIG K0 UNDERFLOWS" , "X", x, '.' ); *bk0 = 0.0; *status = 1; } else { *bk0 = (expf( 10.0 - x )**bk0)*EXPM10; } } } } return; #undef XKN #undef XIN } /* end of function */