slatec-bessel-cpp/f2c/dbsi1e.cpp

/* dbsi1e.f -- translated by f2c (version 20100827).
   This file no longer depends on f2c.
*/

#include "slatec-internal.hpp"

/* Table of constant values */

static integer const c__3 = 3;
static integer const c__17 = 17;
static integer const c__46 = 46;
static integer const c__69 = 69;
static integer const c__1 = 1;

/* Initialized data */

static double const bi1cs[17] = { -.0019717132610998597316138503218149,
                                  .40734887667546480608155393652014,
                                  .034838994299959455866245037783787,
                                  .0015453945563001236038598401058489,
                                  4.188852109837778412945883200412e-5,
                                  7.6490267648362114741959703966069e-7,
                                  1.0042493924741178689179808037238e-8,
                                  9.9322077919238106481371298054863e-11,
                                  7.6638017918447637275200171681349e-13,
                                  4.741418923816739498038809194816e-15,
                                  2.4041144040745181799863172032e-17,
                                  1.0171505007093713649121100799999e-19,
                                  3.6450935657866949458491733333333e-22,
                                  1.1205749502562039344810666666666e-24,
                                  2.9875441934468088832e-27,
                                  6.9732310939194709333333333333333e-30,
                                  1.43679482206208e-32 };
static double const ai1cs[46] = { -.02846744181881478674100372468307,
                                  -.01922953231443220651044448774979,
                                  -6.115185857943788982256249917785e-4,
                                  -2.069971253350227708882823777979e-5,
                                  8.585619145810725565536944673138e-6,
                                  1.04949824671159086251745399786e-6,
                                  -2.918338918447902202093432326697e-7,
                                  -1.559378146631739000160680969077e-8,
                                  1.318012367144944705525302873909e-8,
                                  -1.448423418183078317639134467815e-9,
                                  -2.90851224399314209482504099301e-10,
                                  1.266388917875382387311159690403e-10,
                                  -1.66494777291922067062417839858e-11,
                                  -1.666653644609432976095937154999e-12,
                                  1.242602414290768265232168472017e-12,
                                  -2.731549379672432397251461428633e-13,
                                  2.023947881645803780700262688981e-14,
                                  7.307950018116883636198698126123e-15,
                                  -3.332905634404674943813778617133e-15,
                                  7.17534655851295374354225466567e-16,
                                  -6.982530324796256355850629223656e-17,
                                  -1.299944201562760760060446080587e-17,
                                  8.12094286424279889205467834286e-18,
                                  -2.194016207410736898156266643783e-18,
                                  3.630516170029654848279860932334e-19,
                                  -1.695139772439104166306866790399e-20,
                                  -1.288184829897907807116882538222e-20,
                                  5.694428604967052780109991073109e-21,
                                  -1.459597009090480056545509900287e-21,
                                  2.514546010675717314084691334485e-22,
                                  -1.844758883139124818160400029013e-23,
                                  -6.339760596227948641928609791999e-24,
                                  3.46144110203101111110814662656e-24,
                                  -1.017062335371393547596541023573e-24,
                                  2.149877147090431445962500778666e-25,
                                  -3.045252425238676401746206173866e-26,
                                  5.238082144721285982177634986666e-28,
                                  1.443583107089382446416789503999e-27,
                                  -6.121302074890042733200670719999e-28,
                                  1.700011117467818418349189802666e-28,
                                  -3.596589107984244158535215786666e-29,
                                  5.448178578948418576650513066666e-30,
                                  -2.731831789689084989162564266666e-31,
                                  -1.858905021708600715771903999999e-31,
                                  9.212682974513933441127765333333e-32,
                                  -2.813835155653561106370833066666e-32 };
static double const ai12cs[69] = { .02857623501828012047449845948469,
                                   -.009761097491361468407765164457302,
                                   -1.105889387626237162912569212775e-4,
                                   -3.882564808877690393456544776274e-6,
                                   -2.512236237870208925294520022121e-7,
                                   -2.631468846889519506837052365232e-8,
                                   -3.835380385964237022045006787968e-9,
                                   -5.589743462196583806868112522229e-10,
                                   -1.897495812350541234498925033238e-11,
                                   3.252603583015488238555080679949e-11,
                                   1.412580743661378133163366332846e-11,
                                   2.03562854414708950722452613684e-12,
                                   -7.198551776245908512092589890446e-13,
                                   -4.083551111092197318228499639691e-13,
                                   -2.101541842772664313019845727462e-14,
                                   4.272440016711951354297788336997e-14,
                                   1.042027698412880276417414499948e-14,
                                   -3.814403072437007804767072535396e-15,
                                   -1.880354775510782448512734533963e-15,
                                   3.308202310920928282731903352405e-16,
                                   2.962628997645950139068546542052e-16,
                                   -3.209525921993423958778373532887e-17,
                                   -4.650305368489358325571282818979e-17,
                                   4.414348323071707949946113759641e-18,
                                   7.517296310842104805425458080295e-18,
                                   -9.314178867326883375684847845157e-19,
                                   -1.242193275194890956116784488697e-18,
                                   2.414276719454848469005153902176e-19,
                                   2.026944384053285178971922860692e-19,
                                   -6.394267188269097787043919886811e-20,
                                   -3.049812452373095896084884503571e-20,
                                   1.612841851651480225134622307691e-20,
                                   3.56091396430992505451027090462e-21,
                                   -3.752017947936439079666828003246e-21,
                                   -5.787037427074799345951982310741e-23,
                                   7.759997511648161961982369632092e-22,
                                   -1.452790897202233394064459874085e-22,
                                   -1.318225286739036702121922753374e-22,
                                   6.116654862903070701879991331717e-23,
                                   1.376279762427126427730243383634e-23,
                                   -1.690837689959347884919839382306e-23,
                                   1.430596088595433153987201085385e-24,
                                   3.409557828090594020405367729902e-24,
                                   -1.309457666270760227845738726424e-24,
                                   -3.940706411240257436093521417557e-25,
                                   4.277137426980876580806166797352e-25,
                                   -4.424634830982606881900283123029e-26,
                                   -8.734113196230714972115309788747e-26,
                                   4.045401335683533392143404142428e-26,
                                   7.067100658094689465651607717806e-27,
                                   -1.249463344565105223002864518605e-26,
                                   2.867392244403437032979483391426e-27,
                                   2.04429289250429267028177957421e-27,
                                   -1.518636633820462568371346802911e-27,
                                   8.110181098187575886132279107037e-29,
                                   3.58037935477358609112717370327e-28,
                                   -1.692929018927902509593057175448e-28,
                                   -2.222902499702427639067758527774e-29,
                                   5.424535127145969655048600401128e-29,
                                   -1.787068401578018688764912993304e-29,
                                   -6.56547906872281493882392943788e-30,
                                   7.807013165061145280922067706839e-30,
                                   -1.816595260668979717379333152221e-30,
                                   -1.287704952660084820376875598959e-30,
                                   1.114548172988164547413709273694e-30,
                                   -1.808343145039336939159368876687e-31,
                                   -2.231677718203771952232448228939e-31,
                                   1.619029596080341510617909803614e-31,
                                   -1.83407990880494141390130843921e-32 };

static float const eta = (float) d1mach_(3) * (float).1;
static integer const nti1 = initds_(bi1cs, &c__17, &eta);
static double const xmin = d1mach_(1) * 2.;
static double const xsml = sqrt(d1mach_(3) * 4.5);
static integer const ntai1 = initds_(ai1cs, &c__46, &eta);
static integer const ntai12 = initds_(ai12cs, &c__69, &eta);

double dbsi1e_(double const *x)
{
    /* System generated locals */
    double ret_val, d__1;

    /* Local variables */
    double y;

/* ***BEGIN PROLOGUE  DBSI1E */
/* ***PURPOSE  Compute the exponentially scaled modified (hyperbolic) */
/*            Bessel function of the first kind of order one. */
/* ***LIBRARY   SLATEC (FNLIB) */
/* ***CATEGORY  C10B1 */
/* ***TYPE      DOUBLE PRECISION (BESI1E-S, DBSI1E-D) */
/* ***KEYWORDS  EXPONENTIALLY SCALED, FIRST KIND, FNLIB, */
/*             HYPERBOLIC BESSEL FUNCTION, MODIFIED BESSEL FUNCTION, */
/*             ORDER ONE, SPECIAL FUNCTIONS */
/* ***AUTHOR  Fullerton, W., (LANL) */
/* ***DESCRIPTION */

/* DBSI1E(X) calculates the double precision exponentially scaled */
/* modified (hyperbolic) Bessel function of the first kind of order */
/* one for double precision argument X.  The result is I1(X) */
/* multiplied by EXP(-ABS(X)). */

/* 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 */

/* 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 */

/* ***REFERENCES  (NONE) */
/* ***ROUTINES CALLED  D1MACH, DCSEVL, INITDS, XERMSG */
/* ***REVISION HISTORY  (YYMMDD) */
/*   770701  DATE WRITTEN */
/*   890531  Changed all specific intrinsics to generic.  (WRB) */
/*   890531  REVISION DATE from Version 3.2 */
/*   891214  Prologue converted to Version 4.0 format.  (BAB) */
/*   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ) */
/* ***END PROLOGUE  DBSI1E */
/* ***FIRST EXECUTABLE STATEMENT  DBSI1E */
    y = abs(*x);
    if (y > 3.) {
	goto L20;
    }

    ret_val = 0.;
    if (y == 0.) {
	return ret_val;
    }

    if (y <= xmin) {
	xermsg_("SLATEC", "DBSI1E", "ABS(X) SO SMALL I1 UNDERFLOWS", &c__1, &
		c__1, (ftnlen)6, (ftnlen)6, (ftnlen)29);
    }
    if (y > xmin) {
	ret_val = *x * .5;
    }
    if (y > xsml) {
	d__1 = y * y / 4.5 - 1.;
	ret_val = *x * (dcsevl_(&d__1, bi1cs, &nti1) + .875);
    }
    ret_val = exp(-y) * ret_val;
    return ret_val;

L20:
    if (y <= 8.) {
	d__1 = (48. / y - 11.) / 5.;
	ret_val = (dcsevl_(&d__1, ai1cs, &ntai1) + .375) / sqrt(y);
    }
    if (y > 8.) {
	d__1 = 16. / y - 1.;
	ret_val = (dcsevl_(&d__1, ai12cs, &ntai12) + .375) / sqrt(y);
    }
    ret_val = f2c::d_sign(&ret_val, x);

    return ret_val;
} /* dbsi1e_ */