/*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 "sbi1k1.h" #include #include /* PARAMETER translations */ #define EXP10 (.2202646579480671651695790e5) #define EXPM10 (.4539992976248485153559152e-4) #define LSQ2PI (.9189385332046727417803296e0) #define LSQPI2 (.2257913526447274323630975e0) /* end of PARAMETER translations */ void /*FUNCTION*/ sbi1k1( float x, float *bi1, float *bk1, long want, long *status) { long int _l0; static long int ntai1, ntai12, ntak1, ntak12, ntk1; float i1ns, y, z; static float boundk, ximax, ximin, xisml, xkmax, xkmin, xksml; static float bi1cs[17]={-.19717132610998597316138503218149e-2, .40734887667546480608155393652014e0,.34838994299959455866245037783787e-1, .15453945563001236038598401058489e-2,.41888521098377784129458832004120e-4, .76490267648362114741959703966069e-6,.10042493924741178689179808037238e-7, .99322077919238106481371298054863e-10,.76638017918447637275200171681349e-12, .47414189238167394980388091948160e-14,.24041144040745181799863172032000e-16, .10171505007093713649121100799999e-18,.36450935657866949458491733333333e-21, .11205749502562039344810666666666e-23,.29875441934468088832000000000000e-26, .69732310939194709333333333333333e-29,.14367948220620800000000000000000e-31}; static float ai1cs[46]={-.2846744181881478674100372468307e-1,-.1922953231443220651044448774979e-1, -.6115185857943788982256249917785e-3,-.2069971253350227708882823777979e-4, .8585619145810725565536944673138e-5,.1049498246711590862517453997860e-5, -.2918338918447902202093432326697e-6,-.1559378146631739000160680969077e-7, .1318012367144944705525302873909e-7,-.1448423418183078317639134467815e-8, -.2908512243993142094825040993010e-9,.1266388917875382387311159690403e-9, -.1664947772919220670624178398580e-10,-.1666653644609432976095937154999e-11, .1242602414290768265232168472017e-11,-.2731549379672432397251461428633e-12, .2023947881645803780700262688981e-13,.7307950018116883636198698126123e-14, -.3332905634404674943813778617133e-14,.7175346558512953743542254665670e-15, -.6982530324796256355850629223656e-16,-.1299944201562760760060446080587e-16, .8120942864242798892054678342860e-17,-.2194016207410736898156266643783e-17, .3630516170029654848279860932334e-18,-.1695139772439104166306866790399e-19, -.1288184829897907807116882538222e-19,.5694428604967052780109991073109e-20, -.1459597009090480056545509900287e-20,.2514546010675717314084691334485e-21, -.1844758883139124818160400029013e-22,-.6339760596227948641928609791999e-23, .3461441102031011111108146626560e-23,-.1017062335371393547596541023573e-23, .2149877147090431445962500778666e-24,-.3045252425238676401746206173866e-25, .5238082144721285982177634986666e-27,.1443583107089382446416789503999e-26, -.6121302074890042733200670719999e-27,.1700011117467818418349189802666e-27, -.3596589107984244158535215786666e-28,.5448178578948418576650513066666e-29, -.2731831789689084989162564266666e-30,-.1858905021708600715771903999999e-30, .9212682974513933441127765333333e-31,-.2813835155653561106370833066666e-31}; static float ai12cs[69]={.2857623501828012047449845948469e-1,-.9761097491361468407765164457302e-2, -.1105889387626237162912569212775e-3,-.3882564808877690393456544776274e-5, -.2512236237870208925294520022121e-6,-.2631468846889519506837052365232e-7, -.3835380385964237022045006787968e-8,-.5589743462196583806868112522229e-9, -.1897495812350541234498925033238e-10,.3252603583015488238555080679949e-10, .1412580743661378133163366332846e-10,.2035628544147089507224526136840e-11, -.7198551776245908512092589890446e-12,-.4083551111092197318228499639691e-12, -.2101541842772664313019845727462e-13,.4272440016711951354297788336997e-13, .1042027698412880276417414499948e-13,-.3814403072437007804767072535396e-14, -.1880354775510782448512734533963e-14,.3308202310920928282731903352405e-15, .2962628997645950139068546542052e-15,-.3209525921993423958778373532887e-16, -.4650305368489358325571282818979e-16,.4414348323071707949946113759641e-17, .7517296310842104805425458080295e-17,-.9314178867326883375684847845157e-18, -.1242193275194890956116784488697e-17,.2414276719454848469005153902176e-18, .2026944384053285178971922860692e-18,-.6394267188269097787043919886811e-19, -.3049812452373095896084884503571e-19,.1612841851651480225134622307691e-19, .3560913964309925054510270904620e-20,-.3752017947936439079666828003246e-20, -.5787037427074799345951982310741e-22,.7759997511648161961982369632092e-21, -.1452790897202233394064459874085e-21,-.1318225286739036702121922753374e-21, .6116654862903070701879991331717e-22,.1376279762427126427730243383634e-22, -.1690837689959347884919839382306e-22,.1430596088595433153987201085385e-23, .3409557828090594020405367729902e-23,-.1309457666270760227845738726424e-23, -.3940706411240257436093521417557e-24,.4277137426980876580806166797352e-24, -.4424634830982606881900283123029e-25,-.8734113196230714972115309788747e-25, .4045401335683533392143404142428e-25,.7067100658094689465651607717806e-26, -.1249463344565105223002864518605e-25,.2867392244403437032979483391426e-26, .2044292892504292670281779574210e-26,-.1518636633820462568371346802911e-26, .8110181098187575886132279107037e-28,.3580379354773586091127173703270e-27, -.1692929018927902509593057175448e-27,-.2222902499702427639067758527774e-28, .5424535127145969655048600401128e-28,-.1787068401578018688764912993304e-28, -.6565479068722814938823929437880e-29,.7807013165061145280922067706839e-29, -.1816595260668979717379333152221e-29,-.1287704952660084820376875598959e-29, .1114548172988164547413709273694e-29,-.1808343145039336939159368876687e-30, -.2231677718203771952232448228939e-30,.1619029596080341510617909803614e-30, -.1834079908804941413901308439210e-31}; static float bk1cs[16]={.25300227338947770532531120868533e-1,-.35315596077654487566723831691801e0, -.12261118082265714823479067930042e0,-.69757238596398643501812920296083e-2, -.17302889575130520630176507368979e-3,-.24334061415659682349600735030164e-5, -.22133876307347258558315252545126e-7,-.14114883926335277610958330212608e-9, -.66669016941993290060853751264373e-12,-.24274498505193659339263196864853e-14, -.70238634793862875971783797120000e-17,-.16543275155100994675491029333333e-19, -.32338347459944491991893333333333e-22,-.53312750529265274999466666666666e-25, -.75130407162157226666666666666666e-28,-.91550857176541866666666666666666e-31}; static float ak1cs[38]={.27443134069738829695257666227266e0,.75719899531993678170892378149290e-1, -.14410515564754061229853116175625e-2,.66501169551257479394251385477036e-4, -.43699847095201407660580845089167e-5,.35402774997630526799417139008534e-6, -.33111637792932920208982688245704e-7,.34459775819010534532311499770992e-8, -.38989323474754271048981937492758e-9,.47208197504658356400947449339005e-10, -.60478356628753562345373591562890e-11,.81284948748658747888193837985663e-12, -.11386945747147891428923915951042e-12,.16540358408462282325972948205090e-13, -.24809025677068848221516010440533e-14,.38292378907024096948429227299157e-15, -.60647341040012418187768210377386e-16,.98324256232648616038194004650666e-17, -.16284168738284380035666620115626e-17,.27501536496752623718284120337066e-18, -.47289666463953250924281069568000e-19,.82681500028109932722392050346666e-20, -.14681405136624956337193964885333e-20,.26447639269208245978085894826666e-21, -.48290157564856387897969868800000e-22,.89293020743610130180656332799999e-23, -.16708397168972517176997751466666e-23,.31616456034040694931368618666666e-24, -.60462055312274989106506410666666e-25,.11678798942042732700718421333333e-25, -.22773741582653996232867840000000e-26,.44811097300773675795305813333333e-27, -.88932884769020194062336000000000e-28,.17794680018850275131392000000000e-28, -.35884555967329095821994666666666e-29,.72906290492694257991679999999999e-30, -.14918449845546227073024000000000e-30,.30736573872934276300799999999999e-31}; static float ak12cs[33]={.6379308343739001036600488534102e-1,.2832887813049720935835030284708e-1, -.2475370673905250345414545566732e-3,.5771972451607248820470976625763e-5, -.2068939219536548302745533196552e-6,.9739983441381804180309213097887e-8, -.5585336140380624984688895511129e-9,.3732996634046185240221212854731e-10, -.2825051961023225445135065754928e-11,.2372019002484144173643496955486e-12, -.2176677387991753979268301667938e-13,.2157914161616032453939562689706e-14, -.2290196930718269275991551338154e-15,.2582885729823274961919939565226e-16, -.3076752641268463187621098173440e-17,.3851487721280491597094896844799e-18, -.5044794897641528977117282508800e-19,.6888673850418544237018292223999e-20, -.9775041541950118303002132480000e-21,.1437416218523836461001659733333e-21, -.2185059497344347373499733333333e-22,.3426245621809220631645388800000e-23, -.5531064394246408232501248000000e-24,.9176601505685995403782826666666e-25, -.1562287203618024911448746666666e-25,.2725419375484333132349439999999e-26, -.4865674910074827992378026666666e-27,.8879388552723502587357866666666e-28, -.1654585918039257548936533333333e-28,.3145111321357848674303999999999e-29, -.6092998312193127612416000000000e-30,.1202021939369815834623999999999e-30, -.2412930801459408841386666666666e-31}; static long nti1 = 0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Ai12cs = &ai12cs[0] - 1; float *const Ai1cs = &ai1cs[0] - 1; float *const Ak12cs = &ak12cs[0] - 1; float *const Ak1cs = &ak1cs[0] - 1; float *const Bi1cs = &bi1cs[0] - 1; float *const Bk1cs = &bk1cs[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 SBI1K1 Krogh Moved external statement up for mangle. *>> 1996-04-27 SBI1K1 Krogh Changes to use .C. and C%%. *>> 1995-11-28 SBI1K1 Krogh Changes to simplify conversion to C. *>> 1995-11-03 SBI1K1 Krogh Removed blanks in numbers for C conversion. *>> 1994-10-20 SBI1K1 Krogh Changes to use M77CON *>> 1994-04-19 SBI1K1 CLL Edited to make DP & SP files similar. *>> 1990-09-25 SBI1K1 WV Snyder JPL Use Fullerton codes from CMLIB *--S replaces "?": ?BI1K1, ?CSEVL, ?ERM1, ?INITS * * Compute hyperbolic Bessel functions I1 and K1. * Approximations originally produced by Wayne Fullerton, LASL. * * ***** Formal Arguments *********************************** * * X [in] is the argument at which the functions are to be evaluated. * BI1, BK1 [out] are the function values. * WANT [integer,in] indicates the functions to be computed, and their * scaling: * ABS(WANT) = * 1 means compute I1(X) * 2 means compute K1(X) * 3 means compute both of I1(X) and K1(X) * WANT < 0 means compute EXP(-X)*I1(X) and/or EXP(X)*K1(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 K1(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 I1(X) overflowed, * -3 means X was zero or negative and K1(X) is to be computed, * -4 means X was so small K1(X) overflows. * 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 BI1 = the largest representable number; * if status = -3 or -4 then BK1 = 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 BI1 ON THE INTERVAL 0. TO 9.00000E+00 * WITH WEIGHTED ERROR 1.44E-32 * LOG WEIGHTED ERROR 31.84 * SIGNIFICANT FIGURES REQUIRED 31.45 * DECIMAL PLACES REQUIRED 32.46 * */ /* SERIES FOR AI1 ON THE INTERVAL 1.25000E-01 TO 3.33333E-01 * WITH WEIGHTED ERROR 2.81E-32 * LOG WEIGHTED ERROR 31.55 * SIGNIFICANT FIGURES REQUIRED 29.93 * DECIMAL PLACES REQUIRED 32.38 * *++ Save data by elements if ~.C. */ /* SERIES FOR AI12 ON THE INTERVAL 0. TO 1.25000E-01 * WITH WEIGHTED ERROR 1.83E-32 * LOG WEIGHTED ERROR 31.74 * SIGNIFICANT FIGURES REQUIRED 29.97 * DECIMAL PLACES REQUIRED 32.66 * *++ Save data by elements if ~.C. */ /* SERIES FOR BK1 ON THE INTERVAL 0. TO 4.00000E+00 * WITH WEIGHTED ERROR 9.16E-32 * LOG WEIGHTED ERROR 31.04 * SIGNIFICANT FIGURES REQUIRED 30.61 * DECIMAL PLACES REQUIRED 31.64 * */ /* SERIES FOR AK1 ON THE INTERVAL 1.25000E-01 TO 5.00000E-01 * WITH WEIGHTED ERROR 3.07E-32 * LOG WEIGHTED ERROR 31.51 * SIGNIFICANT FIGURES REQUIRED 30.71 * DECIMAL PLACES REQUIRED 32.30 * *++ Save data by elements if ~.C. */ /* SERIES FOR AK12 ON THE INTERVAL 0. TO 1.25000E-01 * WITH WEIGHTED ERROR 2.41E-32 * LOG WEIGHTED ERROR 31.62 * SIGNIFICANT FIGURES REQUIRED 30.25 * DECIMAL PLACES REQUIRED 32.38 * *++ Save data by elements if ~.C. */ /* ***** Statement Functions ******************************** * */ #define XIN(z) ((float)(ximax + 0.5*logf( (z)/powif(1.0 - 3.0/(8.0*\ (z)),2) ))) #define XKN(z) ((float)(xkmax - 0.5*logf( (z)/powif(1.0 + 3.0/(8.0*\ (z)),2) ))) /* ***** Executable Statements ****************************** * */ if (nti1 == 0) { z = 0.1*FLT_EPSILON/FLT_RADIX; boundk = 1.693 - 6.965e-3*logf( z ); sinits( bi1cs, 17, z, &nti1 ); sinits( ai1cs, 46, z, &ntai1 ); sinits( ai12cs, 69, z, &ntai12 ); sinits( bk1cs, 16, z, &ntk1 ); sinits( ak1cs, 38, z, &ntak1 ); sinits( ak12cs, 33, z, &ntak12 ); ximin = 2.0e0*FLT_MIN; xisml = sqrtf( 80.0e0*z ); ximax = logf( FLT_MAX ) + LSQ2PI; ximax = XIN( XIN( ximax ) ); xksml = sqrtf( 40.0e0*z ); xkmax = -logf( FLT_MIN ); xkmin = expf( fmaxf( -xkmax, -logf( FLT_MAX ) ) + 0.01e0 ); xkmax += LSQPI2; xkmax = XKN( XKN( XKN( xkmax ) ) ); } if (want == 0 || labs( want ) > 3) { ierm1( "SBI1K1", -1, 0, "WANT = 0 OR ABS(WANT) > 3", "WANT" , want, '.' ); *status = -1; return; } *status = 0; /* Compute I1 if requested, or if needed for K1 computation. * */ y = fabsf( x ); if (y <= 3.0e0) { if (labs( want ) != 2 || y <= boundk) { if (y < ximin) { /* if (y .ne. 0.0E0) * * call SERM1('SBI1K1',2,-1, * * 'ABS(X) SO SMALL I1 UNDERFLOWS', 'X',x,'.') */ *bi1 = 0.0; /* status = 2 */ } else { if (y <= xisml) { i1ns = 0.5e0*x; } else { i1ns = x*(0.875e0 + scsevl( y*y/4.5e0 - 1.0e0, bi1cs, nti1 )); } if (want < 0) { *bi1 = expf( -y )*i1ns; } else { *bi1 = i1ns; } } } } else if (labs( want ) != 2) { if (y <= 8.0e0) { *bi1 = scsevl( (48.0e0/y - 11.0e0)/5.0e0, ai1cs, ntai1 ); } else { *bi1 = scsevl( 16.0e0/y - 1.0e0, ai12cs, ntai12 ); } *bi1 = signf( (0.375e0 + *bi1)/sqrtf( y ), x ); if (want > 0) { if (y > ximax) { serm1( "SBI1K1", -2, 0, "ABS(X) SO BIG I1 OVERFLOWS" , "X", x, '.' ); *status = -2; *bi1 = FLT_MAX; /* y > ximax => y > xkmax */ *bk1 = 0.0; if (x > 0.0) return; } else { *bi1 = (expf( y - 10.0e0 )**bi1)*EXP10; } } } /* Compute K1 if requested. * */ if (labs( want ) > 1) { if (x <= 0.0e0) { serm1( "SBI1K1", -3, 0, "X IS ZERO OR NEGATIVE", "X", x, '.' ); *status = -3; *bk1 = FLT_MAX; } else if (x <= boundk) { if (x < xkmin) { serm1( "SBI1K1", -4, 0, "X IS SO SMALL K1 OVERFLOWS" , "X", x, '.' ); *status = -4; *bk1 = FLT_MAX; } else { if (x <= xksml) { y = 0.e0; } else { y = x*x; } *bk1 = logf( 0.5e0*x )*i1ns + (0.75e0 + scsevl( 0.5e0* y - 1.0e0, bk1cs, ntk1 ))/x; if (want < 0) *bk1 *= expf( x ); } } else { y = 16.0e0/x; if (x <= 8.0e0) { *bk1 = scsevl( (y - 5.0e0)/3.0e0, ak1cs, ntak1 ); } else { *bk1 = scsevl( y - 1.0e0, ak12cs, ntak12 ); } *bk1 = (1.25e0 + *bk1)/sqrtf( x ); if (want > 0) { if (x > xkmax) { serm1( "SBI1K1", 1, -1, "X SO BIG K1 UNDERFLOWS" , "X", x, '.' ); *bk1 = 0.0e0; *status = 1; } else { *bk1 = (expf( 10.0e0 - x )**bk1)*EXPM10; } } } } return; #undef XKN #undef XIN } /* end of function */