/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:55 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "dbi0k0.h"
#include <float.h>
#include <stdlib.h>
/* PARAMETER translations */
#define EXP10 (.2202646579480671651695790e5)
#define EXPM10 (.4539992976248485153559152e-4)
#define LSQ2PI (.9189385332046727417803296e0)
#define LSQPI2 (.2257913526447274323630975e0)
/* end of PARAMETER translations */
void /*FUNCTION*/ dbi0k0(
double x,
double *bi0,
double *bk0,
long want,
long *status)
{
long int _l0, i;
static long int ntai0, ntai02, ntak0, ntak02;
double sumf, sumg, sump, sumq, y, z;
static double boundk, ximax, xisml, xkmax, xksml;
static double 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 double 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 double 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 double p[6]={5.8599221412826100000e-04,1.3166052564989571850e-01,
1.1999463724910714109e01,4.6850901201934832188e02,5.9169059852270512312e03,
2.4708152720399552679e03};
static double q[2]={-2.4994418972832303646e02,2.1312714303849120380e04};
static double f[4]={-1.6414452837299064100e00,-2.9601657892958843866e02,
-1.7733784684952985886e04,-4.0320340761145482298e05};
static double g[3]={-2.5064972445877992730e02,2.9865713163054025489e04,
-1.6128136304458193998e06};
static double pp[10]={1.1394980557384778174e02,3.6832589957340267940e03,
3.1075408980684392399e04,1.0577068948034021957e05,1.7398867902565686251e05,
1.5097646353289914539e05,7.1557062783764037541e04,1.8321525870183537725e04,
2.3444738764199315021e03,1.1600249425076035558e02};
static double qq[10]={2.0013443064949242491e02,4.4329628889746408858e03,
3.1474655750295278825e04,9.7418829762268075784e04,1.5144644673520157801e05,
1.2689839587977598727e05,5.8824616785857027752e04,1.4847228371802360957e04,
1.8821890840982713696e03,9.2556599177304839811e01};
static double 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 double 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 */
double *const Ai02cs = &ai02cs[0] - 1;
double *const Ai0cs = &ai0cs[0] - 1;
double *const Ak02cs = &ak02cs[0] - 1;
double *const Ak0cs = &ak0cs[0] - 1;
double *const Bi0cs = &bi0cs[0] - 1;
double *const F = &f[0] - 1;
double *const G = &g[0] - 1;
double *const P = &p[0] - 1;
double *const Pp = &pp[0] - 1;
double *const Q = &q[0] - 1;
double *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 DBI0K0 Krogh Moved external statement up for mangle.
*>> 1996-04-27 DBI0K0 Krogh Changes to use .C. and C%%.
*>> 1995-11-28 DBI0K0 Krogh Changes to simplify conversion to C.
*>> 1994-10-20 DBI0K0 Krogh Changes to use M77CON
*>> 1990-10-18 DBI0K0 WV Snyder JPL Use Cody K0 for small X
*>> 1990-09-25 DBI0K0 WV Snyder JPL Use Fullerton codes from CMLIB
*--D 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.00000D+00
* WITH WEIGHTED ERROR 9.51D-34
* LOG WEIGHTED ERROR 33.02
* SIGNIFICANT FIGURES REQUIRED 33.31
* DECIMAL PLACES REQUIRED 33.65
* */
/* SERIES FOR AI0 ON THE INTERVAL 1.25000D-01 TO 3.33333D-01
* WITH WEIGHTED ERROR 2.74D-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.25000D-01
* WITH WEIGHTED ERROR 1.97D-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.25000D-01 TO 5.00000D-01
* WITH WEIGHTED ERROR 2.85D-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.25000D-01
* WITH WEIGHTED ERROR 2.30D-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) ((double)(ximax + 0.5*log( (z)/powi(1.0 + 1.0/(8.0*\
(z)),2) )))
#define XKN(z) ((double)(xkmax - 0.5*log( (z)/powi(1.0 - 1.0/(8.0*\
(z)),2) )))
/* ***** Executable Statements ******************************
* */
if (nti0 <= 0)
{
z = 0.1*DBL_EPSILON/FLT_RADIX;
boundk = 1.757 - 6.22e-3*log( z );
dinits( bi0cs, 18, z, &nti0 );
dinits( ai0cs, 46, z, &ntai0 );
dinits( ai02cs, 69, z, &ntai02 );
dinits( ak0cs, 38, z, &ntak0 );
dinits( ak02cs, 33, z, &ntak02 );
xisml = sqrt( 80.0e0*z );
ximax = log( DBL_MAX ) + LSQ2PI;
ximax = XIN( XIN( ximax ) );
xksml = sqrt( 10.0*z*P[6]/P[5] );
/* xksml = sqrt (10.0*z*min(p(6)/p(5),q(2)/q(1))) */
xkmax = LSQPI2 - log( DBL_MIN );
xkmax = XKN( XKN( XKN( xkmax ) ) );
}
if (want == 0 || labs( want ) > 3)
{
ierm1( "DBI0K0", -1, 0, "WANT = 0 OR ABS(WANT) > 3", "WANT"
, want, '.' );
*status = -1;
return;
}
*status = 0;
/* Compute I0 if requested.
* */
if (labs( want ) != 2)
{
y = fabs( x );
if (y <= 3.0)
{
if (y <= xisml)
{
*bi0 = 1.0e0;
}
else
{
*bi0 = 2.75e0 + dcsevl( y*y/4.5e0 - 1.0e0, bi0cs,
nti0 );
}
if (want < 0)
*bi0 *= exp( -y );
}
else
{
if (y <= 8.0e0)
{
*bi0 = dcsevl( (48.0e0/y - 11.0e0)/5.0e0, ai0cs, ntai0 );
}
else
{
*bi0 = dcsevl( 16.0e0/y - 1.0e0, ai02cs, ntai02 );
}
*bi0 = (0.375e0 + *bi0)/sqrt( y );
if (want > 0)
{
if (y > ximax)
{
derm1( "DBI0K0", -2, 0, "ABS(X) SO BIG I0 OVERFLOWS"
, "X", x, '.' );
*status = -2;
*bi0 = DBL_MAX;
/* y > ximax => y > xkmax */
*bk0 = 0.0;
if (x > 0.0)
return;
}
else
{
*bi0 = (exp( y - 10.0 )**bi0)*EXP10;
}
}
}
}
/* Compute K0 if requested.
* */
if (labs( want ) > 1)
{
if (x <= 0.0e0)
{
derm1( "DBI0K0", -3, 0, "X IS ZERO OR NEGATIVE", "X",
x, '.' );
*status = -3;
*bk0 = DBL_MAX;
}
else if (x <= 1.0)
{
/* 0.0 < x <= 1.0 */
*bk0 = log( 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 *= exp( 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*sqrt( x ));
if (want > 0)
*bk0 *= exp( -x );
}
else
{
/* boundk < x */
y = 16.0e0/x;
if (x <= 8.0e0)
{
*bk0 = dcsevl( (y - 5.0e0)/3.0e0, ak0cs, ntak0 );
}
else
{
*bk0 = dcsevl( y - 1.0e0, ak02cs, ntak02 );
}
*bk0 = (1.25e0 + *bk0)/sqrt( x );
if (want > 0)
{
if (x > xkmax)
{
derm1( "DBI0K0", 1, -2, "X SO BIG K0 UNDERFLOWS"
, "X", x, '.' );
*bk0 = 0.0;
*status = 1;
}
else
{
*bk0 = (exp( 10.0 - x )**bk0)*EXPM10;
}
}
}
}
return;
#undef XKN
#undef XIN
} /* end of function */