/*
    Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).
    You may use, copy, modify this code for any purpose and 
    without fee. You may distribute this ORIGINAL package.
*/

/*
 $Log: dbessel.c,v $
 Revision 1.4  2009/04/21 19:04:42  shang
 removed dbeskv routines since they are from numeric recipes

 Revision 1.3  2009/04/21 15:16:09  shang
 added  dbeskv_nu() routine

 Revision 1.2  2003/10/07 16:22:25  shang
 added modified bessel functions: dbeskv()

 Revision 1.1  2003/06/14 21:34:29  borland
 Added several Double precision BESSEL functions in dbessel.c.

*/

/* Bessel I_0(x) function in double precision */

#include <mdb.h>

double dbesi0(double x)
{
    int k;
    double w, t, y;
    static double a[65] = {
        8.5246820682016865877e-11, 2.5966600546497407288e-9, 
        7.9689994568640180274e-8, 1.9906710409667748239e-6, 
        4.0312469446528002532e-5, 6.4499871606224265421e-4, 
        0.0079012345761930579108, 0.071111111109207045212, 
        0.444444444444724909, 1.7777777777777532045, 
        4.0000000000000011182, 3.99999999999999998, 
        1.0000000000000000001, 
        1.1520919130377195927e-10, 2.2287613013610985225e-9, 
        8.1903951930694585113e-8, 1.9821560631611544984e-6, 
        4.0335461940910133184e-5, 6.4495330974432203401e-4, 
        0.0079013012611467520626, 0.071111038160875566622, 
        0.44444450319062699316, 1.7777777439146450067, 
        4.0000000132337935071, 3.9999999968569015366, 
        1.0000000003426703174, 
        1.5476870780515238488e-10, 1.2685004214732975355e-9, 
        9.2776861851114223267e-8, 1.9063070109379044378e-6, 
        4.0698004389917945832e-5, 6.4370447244298070713e-4, 
        0.0079044749458444976958, 0.071105052411749363882, 
        0.44445280640924755082, 1.7777694934432109713, 
        4.0000055808824003386, 3.9999977081165740932, 
        1.0000004333949319118, 
        2.0675200625006793075e-10, -6.1689554705125681442e-10, 
        1.2436765915401571654e-7, 1.5830429403520613423e-6, 
        4.2947227560776583326e-5, 6.3249861665073441312e-4, 
        0.0079454472840953930811, 0.070994327785661860575, 
        0.44467219586283000332, 1.7774588182255374745, 
        4.0003038986252717972, 3.9998233869142057195, 
        1.0000472932961288324, 
        2.7475684794982708655e-10, -3.8991472076521332023e-9, 
        1.9730170483976049388e-7, 5.9651531561967674521e-7, 
        5.1992971474748995357e-5, 5.7327338675433770752e-4, 
        0.0082293143836530412024, 0.069990934858728039037, 
        0.44726764292723985087, 1.7726685170014087784, 
        4.0062907863712704432, 3.9952750700487845355, 
        1.0016354346654179322
    };
    static double b[70] = {
        6.7852367144945531383e-8, 4.6266061382821826854e-7, 
        6.9703135812354071774e-6, 7.6637663462953234134e-5, 
        7.9113515222612691636e-4, 0.0073401204731103808981, 
        0.060677114958668837046, 0.43994941411651569622, 
        2.7420017097661750609, 14.289661921740860534, 
        59.820609640320710779, 188.78998681199150629, 
        399.8731367825601118, 427.56411572180478514, 
        1.8042097874891098754e-7, 1.2277164312044637357e-6, 
        1.8484393221474274861e-5, 2.0293995900091309208e-4, 
        0.0020918539850246207459, 0.019375315654033949297, 
        0.15985869016767185908, 1.1565260527420641724, 
        7.1896341224206072113, 37.354773811947484532, 
        155.80993164266268457, 489.5211371158540918, 
        1030.9147225169564806, 1093.5883545113746958, 
        4.8017305613187493564e-7, 3.261317843912380074e-6, 
        4.9073137508166159639e-5, 5.3806506676487583755e-4, 
        0.0055387918291051866561, 0.051223717488786549025, 
        0.42190298621367914765, 3.0463625987357355872, 
        18.895299447327733204, 97.915189029455461554, 
        407.13940115493494659, 1274.3088990480582632, 
        2670.9883037012547506, 2815.7166284662544712, 
        1.2789926338424623394e-6, 8.6718263067604918916e-6, 
        1.3041508821299929489e-4, 0.001428224737372747892, 
        0.014684070635768789378, 0.13561403190404185755, 
        1.1152592585977393953, 8.0387088559465389038, 
        49.761318895895479206, 257.2684232313529138, 
        1066.8543146269566231, 3328.3874581009636362, 
        6948.8586598121634874, 7288.4893398212481055, 
        3.409350368197032893e-6, 2.3079025203103376076e-5, 
        3.4691373283901830239e-4, 0.003794994977222908545, 
        0.038974209677945602145, 0.3594948380414878371, 
        2.9522878893539528226, 21.246564609514287056, 
        131.28727387146173141, 677.38107093296675421, 
        2802.3724744545046518, 8718.5731420798254081, 
        18141.348781638832286, 18948.925349296308859
    };
    static double c[45] = {
        2.5568678676452702768e-15, 3.0393953792305924324e-14, 
        6.3343751991094840009e-13, 1.5041298011833009649e-11, 
        4.4569436918556541414e-10, 1.746393051427167951e-8, 
        1.0059224011079852317e-6, 1.0729838945088577089e-4, 
        0.05150322693642527738, 
        5.2527963991711562216e-15, 7.202118481421005641e-15, 
        7.2561421229904797156e-13, 1.482312146673104251e-11, 
        4.4602670450376245434e-10, 1.7463600061788679671e-8, 
        1.005922609132234756e-6, 1.0729838937545111487e-4, 
        0.051503226936437300716, 
        1.3365917359358069908e-14, -1.2932643065888544835e-13, 
        1.7450199447905602915e-12, 1.0419051209056979788e-11, 
        4.58047881980598326e-10, 1.7442405450073548966e-8, 
        1.0059461453281292278e-6, 1.0729837434500161228e-4, 
        0.051503226940658446941, 
        5.3771611477352308649e-14, -1.1396193006413731702e-12, 
        1.2858641335221653409e-11, -5.9802086004570057703e-11, 
        7.3666894305929510222e-10, 1.6731837150730356448e-8, 
        1.0070831435812128922e-6, 1.0729733111203704813e-4, 
        0.051503227360726294675, 
        3.7819492084858931093e-14, -4.8600496888588034879e-13, 
        1.6898350504817224909e-12, 4.5884624327524255865e-11, 
        1.2521615963377513729e-10, 1.8959658437754727957e-8, 
        1.0020716710561353622e-6, 1.073037119856927559e-4, 
        0.05150322383300230775
    };

    w = fabs(x);
    if (w < 8.5) {
        t = w * w * 0.0625;
        k = 13 * ((int) t);
        y = (((((((((((a[k] * t + a[k + 1]) * t + 
            a[k + 2]) * t + a[k + 3]) * t + a[k + 4]) * t + 
            a[k + 5]) * t + a[k + 6]) * t + a[k + 7]) * t + 
            a[k + 8]) * t + a[k + 9]) * t + a[k + 10]) * t + 
            a[k + 11]) * t + a[k + 12];
    } else if (w < 12.5) {
        k = (int) w;
        t = w - k;
        k = 14 * (k - 8);
        y = ((((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12]) * t + b[k + 13];
    } else {
        t = 60 / w;
        k = 9 * ((int) t);
        y = ((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * sqrt(t) * exp(w);
    }
    return y;
}

/* Bessel I_1(x) function in double precision */

double dbesi1(double x)
{
    int k;
    double w, t, y;
    static double a[60] = {
        1.2787464404046789181e-10, 3.5705860060088241077e-9, 
        9.961153761934733504e-8, 2.2395070088633043177e-6, 
        4.0312466928887462346e-5, 5.6437387840203722356e-4, 
        0.0059259259312934746096, 0.04444444444349900887, 
        0.22222222222232042719, 0.66666666666666139867, 
        1.0000000000000001106, 0.49999999999999999962, 
        1.7281952384448634449e-10, 3.064720455997639013e-9, 
        1.0237662138842827028e-7, 2.2299494417341498163e-6, 
        4.0335364374929326943e-5, 5.6433440269141349899e-4, 
        0.0059259754885893798654, 0.04444439941088039787, 
        0.2222222511283502673, 0.66666665422146063244, 
        1.0000000032274936821, 0.49999999961866867205, 
        2.3216048939948030996e-10, 1.7443372702334489579e-9, 
        1.1596478963485415499e-7, 2.1446755518623035147e-6, 
        4.0697440347437076195e-5, 5.6324394900433192204e-4, 
        0.0059283484996093060678, 0.044440673899150997921, 
        0.2222263801685265786, 0.66666358151576732094, 
        1.0000013834029985337, 0.49999971643129650249, 
        3.1013758938255172562e-10, -8.4813676145611694984e-10, 
        1.5544980187411802596e-7, 1.7811109378708045726e-6, 
        4.2945322199060856985e-5, 5.5344850176852353639e-4, 
        0.0059590327716950614802, 0.044371611097707060659, 
        0.22233578241986401111, 0.6665474730046331531, 
        1.0000756505206705927, 0.49997803664415994554, 
        4.1214758313965020365e-10, -5.361331773534742944e-9, 
        2.4661360807517345161e-7, 6.7144593918926723203e-7, 
        5.1988027944945587571e-5, 5.0165568586065803067e-4, 
        0.0061717530047005289953, 0.043745229577317251404, 
        0.22363147971477747996, 0.6647546913111766024, 
        1.0015686689447547657, 0.49941120439785391891
    };
    static double b[70] = {
        6.6324787943143095845e-8, 4.5125928898466638619e-7, 
        6.7937793134877246623e-6, 7.4580507871505926302e-5, 
        7.6866382927334005919e-4, 0.0071185174803491859307, 
        0.058721838073486424416, 0.42473949281714196041, 
        2.6396965606282079123, 13.710008536637016903, 
        57.158647688180932003, 179.46182892089389037, 
        377.57997362398478619, 399.87313678256009819, 
        1.7652713206027939711e-7, 1.1988179244834708057e-6, 
        1.8037851545747139231e-5, 1.9775785516370314656e-4, 
        0.0020354870702829387283, 0.0188221641910322536, 
        0.15500485219010424263, 1.119010001056057321, 
        6.9391565185406617552, 35.948170579648649345, 
        149.41909525103032616, 467.42979492780642582, 
        979.04227423171290408, 1030.9147225169564443, 
        4.7022299276154507603e-7, 3.1878571710170115972e-6, 
        4.7940153875711448496e-5, 5.2496623508411440227e-4, 
        0.0053968661134780824779, 0.049837081920693776234, 
        0.40979593830387765545, 2.9533186922862948404, 
        18.278176130722516369, 94.47649715018912107, 
        391.66075612645333624, 1221.4182034643210345, 
        2548.6177980961291004, 2670.9883037012546541, 
        1.2535083724002034147e-6, 8.484587142065570825e-6, 
        1.2753227372734042108e-4, 0.0013950105363562648921, 
        0.014325473993765291906, 0.13212452778932829125, 
        1.0849287786885151432, 7.8068089156260172673, 
        48.232254570679165833, 248.80659424902394371, 
        1029.0736929484210803, 3200.5629438795801652, 
        6656.7749162019607914, 6948.8586598121632302, 
        3.3439394490599745013e-6, 2.2600596902211837757e-5, 
        3.3955927589987356838e-4, 0.0037105306061050972474, 
        0.038065263634919156421, 0.35068223415665236079, 
        2.8760027832105027316, 20.665999500843274339, 
        127.47939148516390205, 656.43636874254000885, 
        2709.524283793247992, 8407.1174233600734871, 
        17437.146284159740233, 18141.3487816388316
    };
    static double c[45] = {
        -2.8849790431465382128e-15, -3.5125350943844774657e-14, 
        -7.485086701370741975e-13, -1.8383904048277485153e-11, 
        -5.7303556446977223342e-10, -2.4449502737311496525e-8, 
        -1.6765373351766929724e-6, -3.2189516835265773471e-4, 
        0.051503226936425277377, 
        -5.8674306822281631119e-15, -9.4884898451194085565e-15, 
        -8.503386513660036434e-13, -1.8142997866945285736e-11, 
        -5.7340238386338193949e-10, -2.4449138101742183665e-8, 
        -1.6765375646678855842e-6, -3.2189516826945356325e-4, 
        0.051503226936412017608, 
        -1.4723362506764340882e-14, 1.3945147385179042899e-13, 
        -1.9618041857586930923e-12, -1.3343606394065121821e-11, 
        -5.8649674606973244159e-10, -2.4426060539669553778e-8, 
        -1.6765631828366988006e-6, -3.2189515191449587253e-4, 
        0.051503226931820146445, 
        -5.8203519372580372987e-14, 1.2266326995309845825e-12, 
        -1.3921625844526453237e-11, 6.2228025878281625469e-11, 
        -8.8636681342142794023e-10, -2.3661241616744818608e-8, 
        -1.6777870960740520557e-6, -3.2189402882677074318e-4, 
        0.051503226479551959376, 
        -4.5801527369223291722e-14, 6.7998819697143727209e-13, 
        -4.1624857909290468421e-12, -3.2849009406112440998e-11, 
        -3.247827569043111827e-10, -2.5739209934053714983e-8, 
        -1.6730566573215739195e-6, -3.2190010909008684076e-4, 
        0.05150322986693207715
    };

    w = fabs(x);
    if (w < 8.5) {
        t = w * w * 0.0625;
        k = 12 * ((int) t);
        y = (((((((((((a[k] * t + a[k + 1]) * t + 
            a[k + 2]) * t + a[k + 3]) * t + a[k + 4]) * t + 
            a[k + 5]) * t + a[k + 6]) * t + a[k + 7]) * t + 
            a[k + 8]) * t + a[k + 9]) * t + a[k + 10]) * t + 
            a[k + 11]) * w;
    } else if (w < 12.5) {
        k = (int) w;
        t = w - k;
        k = 14 * (k - 8);
        y = ((((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12]) * t + b[k + 13];
    } else {
        t = 60 / w;
        k = 9 * ((int) t);
        y = ((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * sqrt(t) * exp(w);
    }
    return x < 0 ? -y : y;
}

/* Bessel J_0(x) function in double precision */

double dbesj0(double x)
{
    int k;
    double w, t, y, v, theta;
    static double a[8] = {
        -2.3655394e-12, 4.70889868e-10, 
        -6.78167892231e-8, 6.7816840038636e-6, 
        -4.340277777716935e-4, 0.0156249999999992397, 
        -0.2499999999999999638, 0.9999999999999999997
    };
    static double b[65] = {
        6.26681117e-11, -2.2270614428e-9, 
        6.62981656302e-8, -1.6268486502196e-6, 
        3.21978384111685e-5, -5.00523773331583e-4, 
        0.0059060313537449816, -0.0505265323740109701, 
        0.2936432097610503985, -1.0482565081091638637, 
        1.9181123286040428113, -1.13191994752217001, 
        -0.1965480952704682, 
        4.57457332e-11, -1.5814772025e-9, 
        4.55487446311e-8, -1.0735201286233e-6, 
        2.02015179970014e-5, -2.942392368203808e-4, 
        0.0031801987726150648, -0.0239875209742846362, 
        0.1141447698973777641, -0.2766726722823530233, 
        0.1088620480970941648, 0.5136514645381999197, 
        -0.2100594022073706033, 
        3.31366618e-11, -1.1119090229e-9, 
        3.08823040363e-8, -6.956602653104e-7, 
        1.23499947481762e-5, -1.66295194539618e-4, 
        0.0016048663165678412, -0.0100785479932760966, 
        0.0328996815223415274, -0.0056168761733860688, 
        -0.2341096400274429386, 0.2551729256776404262, 
        0.2288438186148935667, 
        2.38007203e-11, -7.731046439e-10, 
        2.06237001152e-8, -4.412291442285e-7, 
        7.3107766249655e-6, -8.91749801028666e-5, 
        7.34165451384135e-4, -0.0033303085445352071, 
        0.0015425853045205717, 0.0521100583113136379, 
        -0.1334447768979217815, -0.1401330292364750968, 
        0.2685616168804818919, 
        1.6935595e-11, -5.308092192e-10, 
        1.35323005576e-8, -2.726650587978e-7, 
        4.151324014176e-6, -4.43353052220157e-5, 
        2.815740758993879e-4, -4.393235121629007e-4, 
        -0.0067573531105799347, 0.0369141914660130814, 
        0.0081673361942996237, -0.257338128589888186, 
        0.0459580257102978932
    };
    static double c[70] = {
        -3.009451757e-11, -1.4958003844e-10, 
        5.06854544776e-9, 1.863564222012e-8, 
        -6.0304249068078e-7, -1.47686259937403e-6, 
        4.714331342682714e-5, 6.286305481740818e-5, 
        -0.00214137170594124344, -8.9157336676889788e-4, 
        0.04508258728666024989, -0.00490362805828762224, 
        -0.27312196367405374426, 0.04193925184293450356, 
        -7.1245356e-12, -4.1170814825e-10, 
        1.38012624364e-9, 5.704447670683e-8, 
        -1.9026363528842e-7, -5.33925032409729e-6, 
        1.736064885538091e-5, 3.0692619152608375e-4, 
        -9.2598938200644367e-4, -0.00917934265960017663, 
        0.02287952522866389076, 0.10545197546252853195, 
        -0.16126443075752985095, -0.19392874768742235538, 
        2.128344556e-11, -3.1053910272e-10, 
        -3.34979293158e-9, 4.50723289505e-8, 
        3.6437959146427e-7, -4.46421436266678e-6, 
        -2.523429344576552e-5, 2.7519882931758163e-4, 
        9.7185076358599358e-4, -0.00898326746345390692, 
        -0.01665959196063987584, 0.11456933464891967814, 
        0.07885001422733148815, -0.23664819446234712621, 
        3.035295055e-11, 5.486066835e-11, 
        -5.01026824811e-9, -5.0124684786e-9, 
        5.8012340163034e-7, 1.6788922416169e-7, 
        -4.373270270147275e-5, 1.183898532719802e-5, 
        0.00189863342862291449, -0.0011375924956163613, 
        -0.03846797195329871681, 0.02389746880951420335, 
        0.22837862066532347461, -0.06765394811166522844, 
        1.279875977e-11, 3.5925958103e-10, 
        -2.28037105967e-9, -4.852770517176e-8, 
        2.8696428000189e-7, 4.40131125178642e-6, 
        -2.366617753349105e-5, -2.4412456252884129e-4, 
        0.00113028178539430542, 0.0070847051391978908, 
        -0.02526914792327618386, -0.08006137953480093426, 
        0.16548380461475971846, 0.14688405470042110229
    };
    static double d[52] = {
        1.059601355592185731e-14, -2.71150591218550377e-13, 
        8.6514809056201638e-12, -4.6264028554286627e-10, 
        5.0815403835647104e-8, -1.76722552048141208e-5, 
        0.16286750396763997378, 2.949651820598278873e-13, 
        -8.818215611676125741e-12, 3.571119876162253451e-10, 
        -2.63192412099371706e-8, 4.709502795656698909e-6, 
        -0.005208333333333283282, 
        7.18344107717531977e-15, -2.51623725588410308e-13, 
        8.6017784918920604e-12, -4.6256876614290359e-10, 
        5.0815343220437937e-8, -1.7672255176494197e-5, 
        0.16286750396763433767, 2.2327570859680094777e-13, 
        -8.464594853517051292e-12, 3.563766464349055183e-10, 
        -2.631843986737892965e-8, 4.70950234228865941e-6, 
        -0.0052083333332278466225, 
        5.15413392842889366e-15, -2.27740238380640162e-13, 
        8.4827767197609014e-12, -4.6224753682737618e-10, 
        5.0814848128929134e-8, -1.7672254763876748e-5, 
        0.16286750396748926663, 1.7316195320192170887e-13, 
        -7.971122772293919646e-12, 3.544039469911895749e-10, 
        -2.631443902081701081e-8, 4.709498228695400603e-6, 
        -0.005208333331514365361, 
        3.84653681453798517e-15, -2.04464520778789011e-13, 
        8.3089298605177838e-12, -4.6155016158412096e-10, 
        5.081326369646665e-8, -1.76722528311426167e-5, 
        0.1628675039665006593, 1.3797879972460878797e-13, 
        -7.448089381011684812e-12, 3.51273379710695978e-10, 
        -2.630500895563592722e-8, 4.709483934775839193e-6, 
        -0.0052083333227940760113
    };

    w = fabs(x);
    if (w < 1) {
        t = w * w;
        y = ((((((a[0] * t + a[1]) * t + 
            a[2]) * t + a[3]) * t + a[4]) * t + 
            a[5]) * t + a[6]) * t + a[7];
    } else if (w < 8.5) {
        t = w * w * 0.0625;
        k = (int) t;
        t -= k + 0.5;
        k *= 13;
        y = (((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12];
    } else if (w < 12.5) {
        k = (int) w;
        t = w - (k + 0.5);
        k = 14 * (k - 8);
        y = ((((((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * t + c[k + 9]) * t + c[k + 10]) * t + 
            c[k + 11]) * t + c[k + 12]) * t + c[k + 13];
    } else {
        v = 24 / w;
        t = v * v;
        k = 13 * ((int) t);
        y = ((((((d[k] * t + d[k + 1]) * t + 
            d[k + 2]) * t + d[k + 3]) * t + d[k + 4]) * t + 
            d[k + 5]) * t + d[k + 6]) * sqrt(v);
        theta = (((((d[k + 7] * t + d[k + 8]) * t + 
            d[k + 9]) * t + d[k + 10]) * t + d[k + 11]) * t + 
            d[k + 12]) * v - 0.78539816339744830962;
        y *= cos(w + theta);
    }
    return y;
}

/* Bessel J_1(x) function in double precision */

double dbesj1(double x)
{
    int k;
    double w, t, y, v, theta;
    static double a[8] = {
        -1.4810349e-13, 3.363594618e-11, 
        -5.65140051697e-9, 6.7816840144764e-7, 
        -5.425347222188379e-5, 0.00260416666666662438, 
        -0.06249999999999999799, 0.49999999999999999998
    };
    static double b[65] = {
        2.43721316e-12, -9.400554763e-11, 
        3.0605338998e-9, -8.287270492518e-8, 
        1.83020515991344e-6, -3.219783841164382e-5, 
        4.3795830161515318e-4, -0.00442952351530868999, 
        0.03157908273375945955, -0.14682160488052520107, 
        0.39309619054093640008, -0.4795280821510107028, 
        0.1414899934402712514, 
        1.82119257e-12, -6.862117678e-11, 
        2.1732790836e-9, -5.69359291782e-8, 
        1.20771046483277e-6, -2.020151799736374e-5, 
        2.5745933218048448e-4, -0.00238514907946126334, 
        0.01499220060892984289, -0.05707238494868888345, 
        0.10375225210588234727, -0.02721551202427354117, 
        -0.06420643306727498985, 
        1.352611196e-12, -4.9706947875e-11, 
        1.527944986332e-9, -3.8602878823401e-8, 
        7.82618036237845e-7, -1.23499947484511e-5, 
        1.45508295194426686e-4, -0.001203649737425854162, 
        0.006299092495799005109, -0.016449840761170764763, 
        0.002106328565019748701, 0.05852741000686073465, 
        -0.031896615709705053191, 
        9.97982124e-13, -3.5702556073e-11, 
        1.062332772617e-9, -2.5779624221725e-8, 
        4.96382962683556e-7, -7.310776625173004e-6, 
        7.8028107569541842e-5, -5.50624088538081113e-4, 
        0.002081442840335570371, -7.71292652260286633e-4, 
        -0.019541271866742634199, 0.033361194224480445382, 
        0.017516628654559387164, 
        7.31050661e-13, -2.5404499912e-11, 
        7.29360079088e-10, -1.6915375004937e-8, 
        3.06748319652546e-7, -4.151324014331739e-6, 
        3.8793392054271497e-5, -2.11180556924525773e-4, 
        2.74577195102593786e-4, 0.003378676555289966782, 
        -0.013842821799754920148, -0.002041834048574905921, 
        0.032167266073736023299
    };
    static double c[70] = {
        -1.185964494e-11, 3.9110295657e-10, 
        1.80385519493e-9, -5.575391345723e-8, 
        -1.8635897017174e-7, 5.42738239401869e-6, 
        1.181490114244279e-5, -3.300031939852107e-4, 
        -3.7717832892725053e-4, 0.01070685852970608288, 
        0.00356629346707622489, -0.13524776185998074716, 
        0.00980725611657523952, 0.27312196367405374425, 
        -3.029591097e-11, 9.259293559e-11, 
        4.96321971223e-9, -1.518137078639e-8, 
        -5.7045127595547e-7, 1.71237271302072e-6, 
        4.271400348035384e-5, -1.2152454198713258e-4, 
        -0.00184155714921474963, 0.00462994691003219055, 
        0.03671737063840232452, -0.06863857568599167175, 
        -0.21090395092505707655, 0.16126443075752985095, 
        -2.19760208e-11, -2.7659100729e-10, 
        3.74295124827e-9, 3.684765777023e-8, 
        -4.5072801091574e-7, -3.27941630669276e-6, 
        3.5713715545163e-5, 1.7664005411843533e-4, 
        -0.00165119297594774104, -0.00485925381792986774, 
        0.03593306985381680131, 0.04997877588191962563, 
        -0.22913866929783936544, -0.07885001422733148814, 
        5.16292316e-12, -3.9445956763e-10, 
        -6.6220021263e-10, 5.511286218639e-8, 
        5.01257940078e-8, -5.22111059203425e-6, 
        -1.34311394455105e-6, 3.0612891890766805e-4, 
        -7.103391195326182e-5, -0.00949316714311443491, 
        0.00455036998246516948, 0.11540391585989614784, 
        -0.04779493761902840455, -0.2283786206653234746, 
        2.697817493e-11, -1.6633326949e-10, 
        -4.3313486035e-9, 2.508404686362e-8, 
        4.8528284780984e-7, -2.58267851112118e-6, 
        -3.521049080466759e-5, 1.6566324273339952e-4, 
        0.00146474737522491617, -0.00565140892697147306, 
        -0.028338820556793004, 0.07580744376982855057, 
        0.16012275906960187978, -0.16548380461475971845
    };
    static double d[52] = {
        -1.272346002224188092e-14, 3.370464692346669075e-13, 
        -1.144940314335484869e-11, 6.863141561083429745e-10, 
        -9.491933932960924159e-8, 5.301676561445687562e-5, 
        0.162867503967639974, -3.652982212914147794e-13, 
        1.151126750560028914e-11, -5.165585095674343486e-10, 
        4.657991250060549892e-8, -1.186794704692706504e-5, 
        0.01562499999999994026, 
        -8.713069680903981555e-15, 3.140780373478474935e-13, 
        -1.139089186076256597e-11, 6.862299023338785566e-10, 
        -9.491926788274594674e-8, 5.301676558106268323e-5, 
        0.162867503967646622, -2.792555727162752006e-13, 
        1.108650207651756807e-11, -5.156745588549830981e-10, 
        4.657894859077370979e-8, -1.186794650130550256e-5, 
        0.01562499999987299901, 
        -6.304859171204770696e-15, 2.857249044208791652e-13, 
        -1.124956921556753188e-11, 6.858482894906716661e-10, 
        -9.49186795351689846e-8, 5.301676509057781574e-5, 
        0.1628675039678191167, -2.185193490132496053e-13, 
        1.048820673697426074e-11, -5.132819367467680132e-10, 
        4.65740943737299422e-8, -1.186794150862988921e-5, 
        0.01562499999779270706, 
        -4.74041720979200985e-15, 2.578715253644144182e-13, 
        -1.104148898414138857e-11, 6.850134201626289183e-10, 
        -9.49167823417491964e-8, 5.301676277588728159e-5, 
        0.1628675039690033136, -1.75512205749384229e-13, 
        9.848723331445182397e-12, -5.094535425482245697e-10, 
        4.656255982268609304e-8, -1.186792402114394891e-5, 
        0.01562499998712198636
    };

    w = fabs(x);
    if (w < 1) {
        t = w * w;
        y = (((((((a[0] * t + a[1]) * t + 
            a[2]) * t + a[3]) * t + a[4]) * t + 
            a[5]) * t + a[6]) * t + a[7]) * w;
    } else if (w < 8.5) {
        t = w * w * 0.0625;
        k = (int) t;
        t -= k + 0.5;
        k *= 13;
        y = ((((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12]) * w;
    } else if (w < 12.5) {
        k = (int) w;
        t = w - (k + 0.5);
        k = 14 * (k - 8);
        y = ((((((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * t + c[k + 9]) * t + c[k + 10]) * t + 
            c[k + 11]) * t + c[k + 12]) * t + c[k + 13];
    } else {
        v = 24 / w;
        t = v * v;
        k = 13 * ((int) t);
        y = ((((((d[k] * t + d[k + 1]) * t + 
            d[k + 2]) * t + d[k + 3]) * t + d[k + 4]) * t + 
            d[k + 5]) * t + d[k + 6]) * sqrt(v);
        theta = (((((d[k + 7] * t + d[k + 8]) * t + 
            d[k + 9]) * t + d[k + 10]) * t + d[k + 11]) * t + 
            d[k + 12]) * v - 0.78539816339744830962;
        y *= sin(w + theta);
    }
    return x < 0 ? -y : y;
}

/* Bessel K_0(x) function in double precision */

double dbesk0(double x)
{
    int k;
    double t, y;
    static double a[16] = {
        2.4307270476772195953e-12, 4.709166636378530437e-10, 
        6.7816861334344265568e-8, 6.7816840204737508252e-6, 
        4.3402777777915334676e-4, 0.015624999999999872796, 
        0.25000000000000000448, 0.99999999999999999997, 
        6.5878327432224993071e-12, 1.2083308769932888218e-9, 
        1.6271062073716412046e-7, 1.4914719278555277887e-5, 
        8.4603509071212245667e-4, 0.02524892993216233391, 
        0.27898287891460312491, 0.11593151565841244874
    };
    static double b[112] = {
        -4.6430702971053162197e-13, 1.037793605956372823e-11, 
        -1.0298475936392057807e-10, 5.3632747492333959219e-10, 
        -2.1674628861036068105e-10, -2.3316071545820437669e-8, 
        2.2557819578691704059e-7, -9.2325694638587080009e-7, 
        -3.3569097781613661759e-6, 8.7355061305812582974e-5, 
        -6.8021202111645760475e-4, 2.7434654781323362319e-4, 
        0.10031787169953909561, 0.42102443824070833334, 
        4.1447451117883103686e-12, -3.4026589638576604315e-11, 
        9.3398790624638977468e-12, 1.518418175079985263e-9, 
        -1.1364911665083029464e-8, 2.0619457602095915719e-8, 
        3.043101803757224363e-7, -2.974973626447455551e-6, 
        8.0143661611467038568e-6, 8.0937525149549218398e-5, 
        -0.0010356346549612699886, 0.0028534806627578638795, 
        0.097369634474060441807, 0.32175066577856452683, 
        1.117088257074072752e-13, -8.2865909408297066068e-11, 
        9.4656678749191182763e-10, -3.583201984184788338e-9, 
        -9.5017955656904252761e-9, 1.5200595674883329093e-7, 
        -3.866326257135605998e-7, -3.3350340828235103499e-6, 
        2.9359886663960844231e-5, -1.1266401822556801563e-5, 
        -0.0012113572742435576205, 0.0063158973673701376253, 
        0.088291790250128171341, 0.22833982383240512262, 
        -3.2880638807053948433e-11, 4.3194884830465283512e-10, 
        -1.7455089683104033093e-9, -3.2437330799994764516e-9, 
        4.7393655539139519778e-8, -1.1929265603456272466e-8, 
        -1.3177845881013419388e-6, 3.3873375636197969526e-6, 
        3.2729835880668256625e-5, -1.8367283883002494561e-4, 
        -8.2830996454188084408e-4, 0.0095512732229514251931, 
        0.072233832113719266702, 0.14753187103603405298, 
        7.9998492614150860098e-11, -7.025734670268613949e-10, 
        7.889882162708458627e-10, 1.1294796399671507085e-8, 
        -1.1360539648638059137e-8, -3.0346309115270564487e-7, 
        3.2235585426189451721e-7, 8.3575612102298214948e-6, 
        -8.5169628089198208211e-6, -2.5740175232173357342e-4, 
        1.246273401468915277e-4, 0.01068323286919220345, 
        0.051515690033637395779, 0.085465862953544883657, 
        -8.6111506537356531608e-11, 5.1862926131024597823e-10, 
        7.5884324949371110022e-10, -6.4011975813006767417e-9, 
        -4.1966181325111763156e-8, 9.1306285446881485314e-8, 
        1.3573638315827954034e-6, 4.8683213252735694701e-7, 
        -3.8805424608710197066e-5, -1.183898646868898061e-4, 
        9.2796213947750964945e-4, 0.0089611057737319027776, 
        0.031464453915862785606, 0.04426764808753663078, 
        4.4400123834164610288e-11, -1.1411233140911074336e-10, 
        -8.820067070246705983e-10, -1.9686735373323381456e-9, 
        1.9921120728941773855e-8, 1.454397441858483474e-7, 
        1.8238418041265854754e-8, -4.5363700392899066037e-6, 
        -2.1688068222527688542e-5, 4.54960621666870347e-5, 
        0.0010435238076080528284, 0.0058374528996419979931, 
        0.01661121071042545585, 0.020756008367065750538, 
        -6.5166519951106397214e-12, -5.857218285878853958e-11, 
        1.5550375065815375404e-10, 1.9526509484993563229e-9, 
        9.2637123346818426594e-9, -1.4136471501812055943e-8, 
        -4.3024895710889717172e-7, -2.3235612243330592076e-6, 
        4.0380616133862188804e-7, 9.2783767992909743602e-5, 
        7.2964887597817095035e-4, 0.0031316245282223273413, 
        0.0078028233022066428316, 0.0090014807263791058095
    };
    static double c[135] = {
        4.5161032649342790231e-11, -4.2774336988557091369e-11, 
        6.0998467173896677777e-10, 1.9845167242599996944e-9, 
        1.3097678767280215271e-8, 7.4505822268382641286e-8, 
        4.2893920879106814989e-7, 2.3900851955655303104e-6, 
        1.2533473009382380357e-5, 5.9693359063879871983e-5, 
        2.477507066108730458e-4, 8.5106703131389516508e-4, 
        0.0022500105115665788755, 0.00404461344545216346, 
        0.0036910983340425942762, 
        3.5732826433251464989e-12, -3.2906649482312266258e-12, 
        7.0873811190464760555e-11, 2.955132058048417712e-10, 
        2.2776940472505079894e-9, 1.5175463612815010036e-8, 
        9.9462487812170164133e-8, 6.14487577978539011e-7, 
        3.486953188290736075e-6, 1.7615836644757657443e-5, 
        7.6373536037879531886e-5, 2.7098571871205999668e-4, 
        7.3399047381788927036e-4, 0.0013439197177355085297, 
        0.0012439943280131230863, 
        3.6343547836242523646e-13, 9.7997961751276137602e-14, 
        1.0184692699811569047e-11, 6.1495184828957652064e-11, 
        5.0238328349302602543e-10, 3.7498626376004337661e-9, 
        2.6689445483857236307e-8, 1.7591899737346368084e-7, 
        1.0486448307010701679e-6, 5.4986458466257148573e-6, 
        2.4521456351751345323e-5, 8.8900942259143832228e-5, 
        2.448394771406830019e-4, 4.5418248688489693045e-4, 
        4.2479574186923180694e-4, 
        5.2460389348163395857e-14, 7.480206302650350354e-14, 
        2.0012201610651998417e-12, 1.4887306044735163359e-11, 
        1.294670541423294035e-10, 1.0391628915892803144e-9, 
        7.8091180499677328456e-9, 5.3694223626907660084e-8, 
        3.3063914804658509029e-7, 1.7776972424421486506e-6, 
        8.0833148098458320202e-6, 2.975555630444881778e-5, 
        8.2945928349220642178e-5, 1.5536921180500112883e-4, 
        1.4647070522281538711e-4, 
        9.7531436733955514559e-15, 2.4084291220447154982e-14, 
        4.7654956400897494468e-13, 4.0200949504810597783e-12, 
        3.6726577109162191533e-11, 3.0939005665422637601e-10, 
        2.4122848979784500179e-9, 1.7071884462645525505e-8, 
        1.0752238955654933405e-7, 5.8844190041189462347e-7, 
        2.7136083303224014597e-6, 1.0102477728604441135e-5, 
        2.8420490721532571809e-5, 5.3637016379451944413e-5, 
        5.0881312956459247572e-5, 
        2.173204986818937726e-15, 7.2720052142815590531e-15, 
        1.28030837955368201e-13, 1.1696825543787717167e-12, 
        1.1083298191597132094e-11, 9.6536661252658773139e-11, 
        7.7242553835198536397e-10, 5.579836626711057562e-9, 
        3.572134529654341437e-8, 1.9806931547193682466e-7, 
        9.2312964655319555313e-7, 3.4666258590861079959e-6, 
        9.8224698307751177077e-6, 1.8648773453825584428e-5, 
        1.7780062316167651812e-5, 
        5.5012463763851934112e-16, 2.2254763392767319419e-15, 
        3.7187669817701214965e-14, 3.5819585377733489628e-13, 
        3.4866061263191556694e-12, 3.110163345062965291e-11, 
        2.5358235662235617663e-10, 1.8597629779492599046e-9, 
        1.2052654739462999992e-8, 6.7501417351172136833e-8, 
        3.1720052198654584574e-7, 1.1993651363602981832e-6, 
        3.4179130317623363474e-6, 6.5208606745808860158e-6, 
        6.2430205476536771454e-6, 
        1.5225407517829491689e-16, 6.9834820025664405161e-16, 
        1.1380182837138781431e-14, 1.1369488761077196511e-13, 
        1.1291168681618466716e-12, 1.0250757630526871007e-11, 
        8.4765287317253141514e-11, 6.2886627779402596211e-10, 
        4.1142865598366029316e-9, 2.3223773435632014408e-8, 
        1.0985095234166396934e-7, 4.1766260951820336228e-7, 
        1.1958609263543792991e-6, 2.2907574647671878055e-6, 
        2.2008253973114914005e-6, 
        4.4863058691420695911e-17, 2.2437356594371819978e-16, 
        3.6107964803015652759e-15, 3.7031193629853392081e-14, 
        3.7341552790439784371e-13, 3.4355950129497564468e-12, 
        2.8719942600171304499e-11, 2.1499646844509516453e-10, 
        1.4171810843455227171e-9, 8.0501118772875784153e-9, 
        3.8281889106330295876e-8, 1.4621673458431979989e-7, 
        4.2029868696411098586e-7, 8.0785884122023473025e-7, 
        7.7845438614204963209e-7
    };
    static double d[40] = {
        -7.9737703860537066166e-14, 1.9543834380466766627e-12, 
        -4.7230794431646733538e-11, 1.4001773785771252004e-9, 
        -5.4864553020583098585e-8, 3.1601984250143742772e-6, 
        -3.3708783204090252161e-4, 0.16180215937964160437, 
        -5.2593898374798632343e-14, 1.7725913926973236457e-12, 
        -4.6672234858122387294e-11, 1.3991653503828889207e-9, 
        -5.4863400156413929639e-8, 3.1601976099900075541e-6, 
        -3.3708783171335864627e-4, 0.1618021593795843376, 
        -3.6135496189875398132e-14, 1.5466239429618130284e-12, 
        -4.5320259146602122624e-11, 1.3945974109459385552e-9, 
        -5.4853994841172088787e-8, 3.1601858228022739196e-6, 
        -3.370878233999830232e-4, 0.16180215937704286491, 
        -2.5640663123518180635e-14, 1.3288079339404032671e-12, 
        -4.3368537955908371563e-11, 1.3848103653102203186e-9, 
        -5.4824335664256344123e-8, 3.1601315173126153586e-6, 
        -3.370877677903569564e-4, 0.16180215935248373474, 
        -1.8678321325292127767e-14, 1.1354310934105733311e-12, 
        -4.1057197297998608931e-11, 1.369399096129635097e-9, 
        -5.4762428935047089835e-8, 3.1599817092775027963e-6, 
        -3.3708756559715893599e-4, 0.1618021592350814424
    };

    if (x < 0.86) {
        t = x * x;
        y = ((((((a[0] * t + a[1]) * t + 
            a[2]) * t + a[3]) * t + a[4]) * t + 
            a[5]) * t + a[6]) * t + a[7];
        y = ((((((a[8] * t + a[9]) * t + 
            a[10]) * t + a[11]) * t + a[12]) * t + 
            a[13]) * t + a[14]) * t + a[15] - y * log(x);
    } else if (x < 4.15) {
        t = x - 5 / x;
        k = (int) (t + 5);
        t = (k - 4) - t;
        k *= 14;
        y = ((((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12]) * t + b[k + 13];
    } else if (x < 12.5) {
        k = (int) x;
        t = (k + 1) - x;
        k = 15 * (k - 4);
        y = (((((((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * t + c[k + 9]) * t + c[k + 10]) * t + 
            c[k + 11]) * t + c[k + 12]) * t + c[k + 13]) * t + 
            c[k + 14];
    } else {
        t = 60 / x;
        k = 8 * ((int) t);
        y = (((((((d[k] * t + d[k + 1]) * t + 
            d[k + 2]) * t + d[k + 3]) * t + d[k + 4]) * t + 
            d[k + 5]) * t + d[k + 6]) * t + d[k + 7]) * 
            sqrt(t) * exp(-x);
    }
    return y;
}

/* Bessel K_1(x) function in double precision */

double dbesk1(double x)
{
    int k;
    double t, y, v;
    static double a[16] = {
        1.5151605362537935201e-13, 3.363790951353651035e-11, 
        5.6514041131016827202e-9, 6.7816840255069534052e-7, 
        5.4253472222259226487e-5, 0.0026041666666666637057, 
        0.06250000000000000009, 0.5, 
        -8.9790303384748696588e-11, -1.4029047449249185771e-8, 
        -1.5592893752540998113e-6, -1.1253607018469017569e-4, 
        -0.0046421827665011579173, -0.085370719728648409609, 
        -0.3079657578292062966, 1.0000000000000000004
    };
    static double b[120] = {
        -9.4055461896630579928e-12, 3.1307934665844664773e-11, 
        4.2005295001519243251e-10, -4.1636196779679820012e-9, 
        1.4483026181700966164e-8, 1.1661000205428816914e-8, 
        -3.5023996724943046209e-7, 1.4404279316339005012e-6, 
        5.358156415715824208e-7, -3.5249754038612334639e-5, 
        1.7150324075631641453e-4, -4.1276362611239191024e-5, 
        -0.0046943110979636602591, 0.035085369853392357659, 
        0.20063574339907819159, 
        3.3998989888944034586e-11, 7.1558979072937373055e-11, 
        -2.9226856932927698732e-9, 1.4591620256525610213e-8, 
        -6.6141635609854161666e-9, -1.9991101838984472332e-7, 
        5.9185836628873530572e-7, 1.9562880347358085687e-6, 
        -1.5814366450418102764e-5, 7.6791682910944612028e-6, 
        2.8354678948323983936e-4, -0.0010217932669418690641, 
        -0.0032205661865210048433, 0.043497494842354644077, 
        0.16110284302315089935, 
        -2.0933987679665541827e-10, 7.9503322090520447316e-10, 
        3.8000150948242575774e-9, -2.3076136195585571309e-8, 
        -2.3902880302550799653e-8, 3.1914500937804377478e-7, 
        3.2639909831082417694e-7, -5.3166994792995439449e-6, 
        -3.1109524694269240094e-6, 9.2575906966353273247e-5, 
        7.5063709094147644746e-7, -0.0017416491592625765379, 
        0.0012138560335171676007, 0.045879687144659643175, 
        0.11566544716132846709, 
        3.1582384905164908749e-10, -1.9959561818098999516e-9, 
        8.6959328920030927557e-10, 1.1642778282445577109e-8, 
        4.3552264337818440471e-8, -1.5057982160481803238e-7, 
        -1.0101729117980989857e-6, 7.7002002510177612013e-7, 
        1.9580574235590194233e-5, 1.9358461980242834361e-5, 
        -3.3932339942485532728e-4, -9.3416673584325090073e-4, 
        0.0055800080455912847227, 0.038668683062477179235, 
        0.072651643500517000658, 
        -1.1554749629758510059e-10, 8.2270678758893273006e-10, 
        -5.0211156951551538591e-10, -1.4929179050554858361e-9, 
        -2.7107940791526366702e-8, -4.2204764086705349384e-8, 
        3.7253098167927628867e-7, 2.4374697215363361156e-6, 
        1.414194200690976837e-6, -4.8766389019473918231e-5, 
        -2.1681387247526720978e-4, 2.9325729929653405236e-4, 
        0.0064087534504827239815, 0.026054628289709454356, 
        0.040156431128194184336, 
        2.5506555170746221691e-11, -1.3521164018407978152e-10, 
        -8.3281235274106699399e-11, -9.7764849575562351891e-10, 
        3.4661828409940354542e-9, 3.9760633711791357544e-8, 
        1.590290664550452993e-7, -1.4919441249454941275e-7, 
        -5.3779684992094351263e-6, -2.7513862296246223142e-5, 
        -9.7880089725297162007e-6, 7.0787668964515789714e-4, 
        0.0046968199862345387583, 0.014745740181663320127, 
        0.020048622219583455723, 
        -3.4824483072529265585e-12, 1.5157161810563380451e-12, 
        8.5303859696700686144e-12, 3.3455414203743741076e-10, 
        2.0226016353844285376e-9, 5.312815400326633499e-9, 
        -3.0799322316418042137e-8, -4.4455408272954712128e-7, 
        -2.4293274626893384034e-6, -3.2129079340119038354e-6, 
        5.922540368307538885e-5, 5.6822962576781683532e-4, 
        0.0027152446516406682732, 0.0074075873691848838485, 
        0.0093044450815739269849, 
        -2.7683216166276377232e-13, 3.1986676777610155465e-12, 
        9.4142986954031445666e-12, 6.7934609179456399334e-11, 
        3.485152941147002933e-11, -2.5785248508896551557e-9, 
        -2.8310220027112571258e-8, -1.6384131113072271115e-7, 
        -3.2521663350596379097e-7, 4.038138875762230716e-6, 
        5.1917606978077281001e-5, 3.3420027947470126154e-4, 
        0.0013699550623118247094, 0.0034405619148342271096, 
        0.0041042919106665762794
    };
    static double c[120] = {
        4.5281968025889407937e-12, 1.0806749918195271176e-11, 
        9.6200972728717669027e-11, 5.721422706362526365e-10, 
        3.6077804282954825099e-9, 2.2465236858536681852e-8, 
        1.367696126430873523e-7, 7.9561767489531997361e-7, 
        4.3014380065615550573e-6, 2.092171390555028559e-5, 
        8.8079183950590176926e-5, 3.0549414408830252064e-4, 
        8.1295715613927890473e-4, 0.0014679809476357079195, 
        0.0013439197177355090057, 
        7.6019964430402432637e-13, -2.261619859915827119e-13, 
        1.7904450823779000744e-11, 9.1467054855312232717e-11, 
        7.1378582044879519122e-10, 4.9925255415445769102e-9, 
        3.3767315471315546644e-8, 2.1350774539167751457e-7, 
        1.2314353082655232903e-6, 6.2918685053670619181e-6, 
        2.7493229298777000013e-5, 9.8085825401369821771e-5, 
        2.6670282677770444935e-4, 4.8967895428135985381e-4, 
        4.5418248688489697144e-4, 
        9.4180115230375147213e-14, 7.5943117003734061145e-14, 
        3.0335730243874287654e-12, 2.0202796115462268051e-11, 
        1.6839020189186971198e-10, 1.2907875663127201526e-9, 
        9.354767612586579892e-9, 6.2471974110281880722e-8, 
        3.7585985422997380441e-7, 1.9838348288114906484e-6, 
        8.8884862203671982034e-6, 3.2333259238682810218e-5, 
        8.9266668913380400243e-5, 1.6589185669844051903e-4, 
        1.5536921180500113394e-4, 
        1.5425475332301107271e-14, 2.8674534590132490434e-14, 
        6.5078462279160216936e-13, 5.0939757793961391211e-12, 
        4.497983746074897552e-11, 3.6662925847520171711e-10, 
        2.7848878755089582413e-9, 1.929812005933947782e-8, 
        1.1950323861976892013e-7, 6.4513432758147478287e-7, 
        2.9422095033982461936e-6, 1.0854433321174584937e-5, 
        3.0307433185818899481e-5, 5.684098144306501785e-5, 
        5.3637016379451945253e-5, 
        3.1077953698439839352e-15, 8.6899496170729520378e-15, 
        1.6258562067326054104e-13, 1.4104842571366761537e-12, 
        1.3019455544084110747e-11, 1.1070466372863950239e-10, 
        8.6890603844230597917e-10, 6.1793722175049967488e-9, 
        3.9058865943755615801e-8, 2.1432806981070368523e-7, 
        9.9034657762983230155e-7, 3.6925185861895664251e-6, 
        1.0399877577259449786e-5, 1.9644939661550210015e-5, 
        1.8648773453825584597e-5, 
        7.2831555285336286457e-16, 2.6077534095895783532e-15, 
        4.4881202059263153495e-14, 4.1674329383944385626e-13, 
        3.9760100480223728037e-12, 3.483597635535118301e-11, 
        2.79932542127702497e-10, 2.0286513276830758107e-9, 
        1.3018343087118439152e-8, 7.2315927974997999365e-8, 
        3.3750708681924201599e-7, 1.2688020879407355571e-6, 
        3.5980954090811587848e-6, 6.8358260635246667316e-6, 
        6.5208606745808860557e-6, 
        1.9026412343503745875e-16, 8.0073765508732553766e-16, 
        1.3245754278055523992e-14, 1.2885201653055058502e-13, 
        1.2600129301230402587e-12, 1.1283306843147549277e-11, 
        9.2261481309646814329e-11, 6.7812033168299846818e-10, 
        4.4020645304595102132e-9, 2.4685719238301517679e-8, 
        1.1611886719473112509e-7, 4.3940380936523135466e-7, 
        1.2529878285546791905e-6, 2.3917218527087570384e-6, 
        2.290757464767187816e-6, 
        5.3709522135744366512e-17, 2.5239221050372845433e-16, 
        4.093314714589908336e-15, 4.1152784247617592367e-14, 
        4.0998840572769381012e-13, 3.7319354625807158852e-12, 
        3.0921671702920868014e-11, 2.2975898538634445343e-10, 
        1.5049754445782364328e-9, 8.5030864719789148982e-9, 
        4.025055939111842381e-8, 1.5312755642491878591e-7, 
        4.3865020375297892208e-7, 8.4059737392822153101e-7, 
        8.0785884122023473319e-7
    };
    static double d[40] = {
        9.2371554649979581914e-14, -2.3111336195788410887e-12, 
        5.7728710326649832559e-11, -1.8002298130091372598e-9, 
        7.6810375010517145638e-8, -5.2669973752193823306e-6, 
        0.0010112634961227401357, 0.16180215937964160466, 
        6.1381146507252683381e-14, -2.1034499679806301862e-12, 
        5.7090233460448415278e-11, -1.7990724350642330817e-9, 
        7.6809056078388019946e-8, -5.2669964425290062357e-6, 
        0.001011263495747828339, 0.16180215937970716383, 
        4.2458150578401296419e-14, -1.8435733128339016981e-12, 
        5.5534955081564656595e-11, -1.7938162188526358466e-9, 
        7.6798230945934117807e-8, -5.2669828728791921259e-6, 
        0.0010112634861753356559, 0.16180215938263409582, 
        3.0314798962267007518e-14, -1.5915009905364214455e-12, 
        5.3275907427402047438e-11, -1.7824862013841369751e-9, 
        7.676389035644707581e-8, -5.2669199860465945909e-6, 
        0.0010112634217687349189, 0.16180215941108227283, 
        2.2211515002229271212e-14, -1.3664088221521734796e-12, 
        5.0585177270502341602e-11, -1.7645432205894533462e-9, 
        7.6691805594577373698e-8, -5.2667455286976269634e-6, 
        0.001011263186281097458, 0.16180215954783127877
    };

    if (x < 0.8) {
        t = x * x;
        y = (((((((a[0] * t + a[1]) * t + 
            a[2]) * t + a[3]) * t + a[4]) * t + 
            a[5]) * t + a[6]) * t + a[7]) * x;
        y = (((((((a[8] * t + a[9]) * t + 
            a[10]) * t + a[11]) * t + a[12]) * t + 
            a[13]) * t + a[14]) * t + a[15]) / x + 
            y * log(x);
    } else if (x < 5.5) {
        v = 3 / x;
        t = x - v;
        k = (int) (t + 3);
        t = (k - 2) - t;
        k *= 15;
        y = ((((((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12]) * t + b[k + 13]) * t + 
            b[k + 14]) * v;
    } else if (x < 12.5) {
        k = (int) x;
        t = (k + 1) - x;
        k = 15 * (k - 5);
        y = (((((((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * t + c[k + 9]) * t + c[k + 10]) * t + 
            c[k + 11]) * t + c[k + 12]) * t + c[k + 13]) * t + 
            c[k + 14];
    } else {
        t = 60 / x;
        k = 8 * ((int) t);
        y = (((((((d[k] * t + d[k + 1]) * t + 
            d[k + 2]) * t + d[k + 3]) * t + d[k + 4]) * t + 
            d[k + 5]) * t + d[k + 6]) * t + d[k + 7]) * 
            sqrt(t) * exp(-x);
    }
    return y;
}

/* Bessel Y_0(x) function in double precision */

double dbesy0(double x)
{
    int k;
    double t, y, v, theta;
    static double a[16] = {
        -1.51249795e-12, 2.9979612902e-10, 
        -4.317352912436e-8, 4.31735413787068e-6, 
        -2.763106650893309e-4, 0.0099471839432433894, 
        -0.15915494309189533339, 0.63661977236758134306, 
        4.09490035e-12, -7.6925095943e-10, 
        1.0358472550303e-7, -9.49500519343105e-6, 
        5.3860266685948738e-4, -0.01607396802593822992, 
        0.17760601686906713536, -0.07380429510868722524
    };
    static double b[112] = {
        -2.958527319e-11, -4.799424787e-11, 
        7.0761365903e-10, 7.86916310954e-9, 
        4.40632174633e-8, 1.047420705431e-7, 
        -7.0359581447993e-7, -1.068225180236166e-5, 
        -7.636033201484949e-5, -3.4318294473105478e-4, 
        -5.0019381740765567e-4, 0.00898034680285547482, 
        0.14783488851798565262, 0.01218096458136948754, 
        -3.790394288e-11, -8.3012324083e-10, 
        -4.16344143065e-9, -1.59750317942e-9, 
        1.1897479804282e-7, 8.9462034665663e-7, 
        2.7927785508031e-6, -7.10267818373816e-6, 
        -1.4158631248747538e-4, -8.9512296835766028e-4, 
        -0.00286419212921734238, 0.00448696141199596815, 
        0.16247276853497028694, 0.16806523426268420517, 
        4.9508820408e-10, 2.40256475451e-9, 
        -3.4128702476e-9, -7.926555593644e-8, 
        -2.7510687602104e-7, 8.2850014912691e-7, 
        1.202838962100801e-5, 4.446440998573947e-5, 
        -6.345799532264893e-5, -0.00153636349657759957, 
        -0.00784582749467635199, -0.01091113270659651752, 
        0.15855078058220179592, 0.33112075987570150227, 
        -1.21381163034e-9, -1.85252445637e-9, 
        2.782181715752e-8, 9.983716245339e-8, 
        -5.1382328105461e-7, -4.13045112612025e-6, 
        6.9219807212556e-7, 1.1346881914173059e-4, 
        4.485714262231131e-4, -7.5862241142011713e-4, 
        -0.01330170710405210993, -0.04337207915110066583, 
        0.10708263075506099964, 0.46937172516119821029, 
        1.23954313286e-9, -1.3150647347e-9, 
        -3.208158247443e-8, 6.936208476862e-8, 
        9.6326321605564e-7, -1.86813677726279e-6, 
        -3.296596487362508e-5, -1.027532398838717e-5, 
        8.8296805857821394e-4, 0.00288502916411875494, 
        -0.00981551522729451135, -0.08175499367255420548, 
        -0.0197090243242566078, 0.51958016003260732459, 
        -5.0293922567e-10, 2.5905819115e-9, 
        3.80069015768e-9, -1.257928760583e-7, 
        2.9666593806065e-7, 5.32199410651341e-6, 
        -1.466013674557952e-5, -2.1196179800800308e-4, 
        1.488907515797525e-4, 0.00598670133725335804, 
        0.00917434640154761191, -0.08592436550224288681, 
        -0.1970089491453584694, 0.41202451505149651586, 
        2.689416111e-11, -6.0225798842e-10, 
        7.16632287642e-9, -3.038834206093e-8, 
        -4.9829955736697e-7, 3.6467866107985e-6, 
        2.599276311771812e-5, -1.6374249229498868e-4, 
        -0.00112385259950125883, 0.00342431039748960465, 
        0.03016663018029814938, -0.02432512430353342927, 
        -0.31797371916576712163, 0.14418001571733803493, 
        1.9463221e-12, -1.3451327719e-10, 
        3.36581427068e-9, 2.163920016772e-8, 
        -5.0868502638387e-7, -1.26562372408989e-6, 
        3.554054033687039e-5, 7.490048474688881e-5, 
        -0.0014238004200886133, -0.00355314274445060456, 
        0.03042001885994619102, 0.07365497405664882878, 
        -0.26882214651298246524, -0.16578636480856581212
    };
    static double c[126] = {
        1.4133441141e-10, -1.16239182574e-9, 
        5.60813675548e-9, 4.361111919906e-8, 
        -1.9980743680592e-7, -6.23828267379045e-6, 
        2.659779111113435e-5, 3.4429768949740887e-4, 
        -0.00128058451746064875, -0.01165906901727504569, 
        0.03800023545011044967, 0.13079665132249175616, 
        -0.30099732306965462304, -0.19470500862950453349, 
        2.985070388e-11, -4.3935370028e-10, 
        -2.70094988368e-9, 5.27013936718e-8, 
        3.5596645009533e-7, -5.5939305510166e-6, 
        -2.415386269153308e-5, 3.4989509210648774e-4, 
        9.8204687587121096e-4, -0.01241793885175385487, 
        -0.0139851984097438996, 0.16758045653582919005, 
        0.02375823895638961834, -0.33948059288191103829, 
        3.440543098e-11, -1.492311224e-11, 
        -5.50307983816e-9, 2.88861169685e-9, 
        6.5975576340131e-7, -6.4253433004282e-7, 
        -5.095354686815779e-5, 6.35545499385876e-5, 
        0.00231750171436134067, -0.00344140664586595266, 
        -0.04796153689875432703, 0.06553727330877767204, 
        0.27409127395927545297, -0.17324243491898233567, 
        1.869953108e-11, 3.5678205287e-10, 
        -3.23369443017e-9, -4.932837601739e-8, 
        4.0636932259569e-7, 4.55294558480472e-6, 
        -3.376596464782138e-5, -2.5758322585307069e-4, 
        0.00167416927479443727, 0.00735300543429220468, 
        -0.03904554680817849396, -0.07593187710651205994, 
        0.25912851048611625179, 0.11731328614820863082, 
        -1.168499781e-11, 4.049157062e-10, 
        1.76042185153e-9, -5.796810963324e-8, 
        -1.7664264701199e-7, 5.65232278687159e-6, 
        1.060939366068237e-5, -3.4382320897779425e-4, 
        -2.8786871040915476e-4, 0.0110369660250104046, 
        0.00105742483935856071, -0.13664188676516064192, 
        0.02616867939853747003, 0.27020510536578747599, 
        -3.09714065e-11, 1.0275535029e-10, 
        5.07546062007e-9, -1.697584561546e-8, 
        -5.8222618994516e-7, 1.92483899143098e-6, 
        4.3388818054202e-5, -1.3714074807064326e-4, 
        -0.00184722933402281351, 0.0051736071998208299, 
        0.03611657815299529568, -0.07491163418624572802, 
        -0.20317989938720766824, 0.17121062620272384486, 
        -2.289344046e-11, -2.7762757821e-10, 
        3.9133021423e-9, 3.673304032817e-8, 
        -4.7258696909223e-7, -3.23109113206792e-6, 
        3.749284946160478e-5, 1.7039067180362504e-4, 
        -0.00172637004639815403, -0.00454152531086626043, 
        0.03717220530377418124, 0.04489395902785574534, 
        -0.23370422835726857839, -0.06753037249787639679, 
        4.73299038e-12, -4.0489249115e-10, 
        -5.6747261733e-10, 5.656135976496e-8, 
        3.564150059329e-8, -5.34821220288168e-6, 
        7.581385461348e-8, 3.1190256463318867e-4, 
        -1.4634564376714538e-4, -0.0095821355155047985, 
        0.00624681472076521329, 0.11513529702572439703, 
        -0.05794254714300082167, -0.22523211169118786537, 
        2.732444065e-11, -1.7726475495e-10, 
        -4.37562701999e-9, 2.68153688556e-8, 
        4.8799100503315e-7, -2.76490187540644e-6, 
        -3.513357925824063e-5, 1.7693155138571443e-4, 
        0.00144533233498513085, -0.00599106092630544975, 
        -0.02759437856689911694, 0.07945362316083459396, 
        0.15383825653750118008, -0.17121430684466928734
    };
    static double d[52] = {
        1.059601355592185731e-14, -2.71150591218550377e-13, 
        8.6514809056201638e-12, -4.6264028554286627e-10, 
        5.0815403835647104e-8, -1.76722552048141208e-5, 
        0.16286750396763997378, 2.949651820598278873e-13, 
        -8.818215611676125741e-12, 3.571119876162253451e-10, 
        -2.63192412099371706e-8, 4.709502795656698909e-6, 
        -0.005208333333333283282, 
        7.18344107717531977e-15, -2.51623725588410308e-13, 
        8.6017784918920604e-12, -4.6256876614290359e-10, 
        5.0815343220437937e-8, -1.7672255176494197e-5, 
        0.16286750396763433767, 2.2327570859680094777e-13, 
        -8.464594853517051292e-12, 3.563766464349055183e-10, 
        -2.631843986737892965e-8, 4.70950234228865941e-6, 
        -0.0052083333332278466225, 
        5.15413392842889366e-15, -2.27740238380640162e-13, 
        8.4827767197609014e-12, -4.6224753682737618e-10, 
        5.0814848128929134e-8, -1.7672254763876748e-5, 
        0.16286750396748926663, 1.7316195320192170887e-13, 
        -7.971122772293919646e-12, 3.544039469911895749e-10, 
        -2.631443902081701081e-8, 4.709498228695400603e-6, 
        -0.005208333331514365361, 
        3.84653681453798517e-15, -2.04464520778789011e-13, 
        8.3089298605177838e-12, -4.6155016158412096e-10, 
        5.081326369646665e-8, -1.76722528311426167e-5, 
        0.1628675039665006593, 1.3797879972460878797e-13, 
        -7.448089381011684812e-12, 3.51273379710695978e-10, 
        -2.630500895563592722e-8, 4.709483934775839193e-6, 
        -0.0052083333227940760113
    };

    if (x < 0.85) {
        t = x * x;
        y = ((((((a[0] * t + a[1]) * t + 
            a[2]) * t + a[3]) * t + a[4]) * t + 
            a[5]) * t + a[6]) * t + a[7];
        y = ((((((a[8] * t + a[9]) * t + 
            a[10]) * t + a[11]) * t + a[12]) * t + 
            a[13]) * t + a[14]) * t + a[15] + y * log(x);
    } else if (x < 4.5) {
        t = x - 4 / x;
        k = (int) (t + 4);
        t -= k - 3.5;
        k *= 14;
        y = ((((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12]) * t + b[k + 13];
    } else if (x < 12.5) {
        k = (int) x;
        t = x - (k + 0.5);
        k = 14 * (k - 4);
        y = ((((((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * t + c[k + 9]) * t + c[k + 10]) * t + 
            c[k + 11]) * t + c[k + 12]) * t + c[k + 13];
    } else {
        v = 24 / x;
        t = v * v;
        k = 13 * ((int) t);
        y = ((((((d[k] * t + d[k + 1]) * t + 
            d[k + 2]) * t + d[k + 3]) * t + d[k + 4]) * t + 
            d[k + 5]) * t + d[k + 6]) * sqrt(v);
        theta = (((((d[k + 7] * t + d[k + 8]) * t + 
            d[k + 9]) * t + d[k + 10]) * t + d[k + 11]) * t + 
            d[k + 12]) * v - 0.78539816339744830962;
        y *= sin(x + theta);
    }
    return y;
}

/* Bessel Y_1(x) function in double precision */

double dbesy1(double x)
{
    int k;
    double t, y, v, theta;
    static double a[16] = {
        -9.46499e-14, 2.141432796e-11, 
        -3.5977944347e-9, 4.3173541397185e-7, 
        -3.453883313621946e-5, 0.00165786399054057257, 
        -0.03978873577297383381, 0.31830988618379067154, 
        -5.571931463e-11, 8.93098810365e-9, 
        -9.9267351748541e-7, 7.164268731546281e-5, 
        -0.00295530533604595764, 0.05434868816050718549, 
        -0.1960570906462388432, -0.63661977236758134367
    };
    static double b[112] = {
        2.553685555e-11, 7.384080449e-11, 
        -2.9605446927e-10, -5.2081739106e-9, 
        -3.676554001376e-8, -1.532652533627e-7, 
        -1.3967947147386e-7, 3.88730181348638e-6, 
        4.052240780599446e-5, 2.5133005806772436e-4, 
        0.00109380629039004777, 0.00298157276348996763, 
        4.7188983238582062e-4, -0.19643830264444486066, 
        1.056916126e-11, 6.2160122556e-10, 
        3.85741197273e-9, 7.491469406e-9, 
        -6.306874425799e-8, -6.869833657799e-7, 
        -3.29785491054157e-6, -5.57893716372629e-6, 
        4.634569030169549e-5, 4.9086790937752985e-4, 
        0.00256488729380268318, 0.00822571374211646414, 
        0.01094512145046992425, -0.19159562910398869317, 
        -4.0211480479e-10, -2.29721816531e-9, 
        6.2931855684e-10, 6.509879184236e-8, 
        3.1589903600713e-7, -3.093448289001e-8, 
        -8.45532004628483e-6, -4.930719004112743e-5, 
        -9.865914153735679e-5, 4.7059803263198826e-4, 
        0.00472271367216891688, 0.01915518502418162917, 
        0.0372239136138415874, -0.1693323078896687811, 
        1.10728466147e-9, 2.05608441891e-9, 
        -2.536199750468e-8, -1.1751665325851e-7, 
        3.404534385077e-7, 4.36083390507752e-6, 
        8.34725302600388e-6, -7.154971448642437e-5, 
        -5.1956165332380432e-4, -0.00100811886858456098, 
        0.00436854257164588225, 0.03425644258611852849, 
        0.0907397117201493499, -0.10791597118914788657, 
        -1.17749768422e-9, 1.68403859797e-9, 
        3.471030171126e-8, -6.227816248757e-8, 
        -1.18405395385178e-6, 1.4478514453927e-7, 
        3.630744090017806e-5, 1.0481065165516454e-4, 
        -5.2483993430622472e-4, -0.00407014291708523893, 
        -0.00574311514714191566, 0.03534768318503458091, 
        0.16539360869993020163, 0.01986240425870529979, 
        4.008050552e-10, -3.36044652853e-9, 
        -9.2387130705e-10, 1.7840823064933e-7, 
        -2.2643988836534e-7, -8.08896250262948e-6, 
        -1.69403516892975e-6, 2.6617475073041589e-4, 
        7.2562705080786576e-4, -0.00399208199997748377, 
        -0.02399816920965700841, -0.0092557474455754244, 
        0.2008281402282066699, 0.21040564991987672576, 
        6.692975871e-11, 1.672202705e-10, 
        -1.433097981251e-8, 5.455077276306e-8, 
        1.04954046999028e-6, -3.39572219821681e-6, 
        -5.541577161092338e-5, 4.344113487321887e-5, 
        0.00184798594297387512, 0.0030104150846815045, 
        -0.02789739141907399514, -0.09421628519129373841, 
        0.09950212445529771576, 0.3749697583874640408, 
        4.528583204e-11, 6.6861092459e-10, 
        -9.2464816857e-9, -8.002058622489e-8, 
        8.7578676979688e-7, 6.28456150216647e-6, 
        -4.276239654834506e-5, -3.4630842915911896e-4, 
        8.9431202169021474e-4, 0.01086380107342193819, 
        0.00145563493724746238, -0.14193383817937981478, 
        -0.15148562296986494725, 0.35720232689066512346
    };
    static double c[126] = {
        3.201620498e-10, -1.77431944127e-9, 
        1.367520951871e-8, -6.171755212447e-8, 
        -4.3602779932091e-7, 1.79827173442894e-6, 
        4.990624965075363e-5, -1.8618453816525107e-4, 
        -0.00206578613620093339, 0.00640292258731756063, 
        0.04663627606907821596, -0.11400070635033154332, 
        -0.26159330264498333979, 0.30099732306965462346, 
        -5.9774708e-13, -3.8827648808e-10, 
        5.27269271362e-9, 2.971059814012e-8, 
        -5.2701406514027e-7, -3.20369808906187e-6, 
        4.475144442564492e-5, 1.6907703884518866e-4, 
        -0.00209937055264007564, -0.00491023437935628693, 
        0.0496717554070154514, 0.04195559522923170315, 
        -0.33516091307165838036, -0.02375823895638961835, 
        2.83380103e-12, -4.4715450235e-10, 
        1.769519961e-10, 6.05337984007e-8, 
        -2.888550814408e-8, -5.93780185020349e-6, 
        5.14027455732102e-6, 3.5667482807472357e-4, 
        -3.813272996260773e-4, -0.01158750857180657934, 
        0.01376562658346365929, 0.14388461069626297876, 
        -0.1310745466175553429, -0.27409127395927545296, 
        2.683697058e-11, -2.4302161327e-10, 
        -4.30151236236e-9, 3.557058903199e-8, 
        4.9328952592926e-7, -3.65732389065382e-6, 
        -3.642356546467712e-5, 2.3636175253326716e-4, 
        0.00154549935517002108, -0.00837084637397210913, 
        -0.02941202173717025197, 0.11713664042453548043, 
        0.15186375421302413107, -0.25912851048611625179, 
        2.914434604e-11, 1.5185121893e-10, 
        -4.88084673389e-9, -1.936460341189e-8, 
        5.7968735781303e-7, 1.58978381365925e-6, 
        -4.521858314881102e-5, -7.426575562367423e-5, 
        0.00206293925392279862, 0.00143934355204571641, 
        -0.04414786410004317489, -0.00317227451807568106, 
        0.273283773530321296, -0.02616867939853747003, 
        6.44550655e-12, 4.0249587279e-10, 
        -1.23789833341e-9, -5.58299757877e-8, 
        1.6975984093144e-7, 5.24003568623094e-6, 
        -1.539871212028104e-5, -3.0372172637669852e-4, 
        8.2284448843625174e-4, 0.00923614667011392611, 
        -0.02069442879928366384, -0.10834973445898588438, 
        0.14982326837249145873, 0.20317989938720766823, 
        -2.115691215e-11, 2.9752035314e-10, 
        3.34739862269e-9, -4.304625868398e-8, 
        -3.6733494871207e-7, 4.25328270524107e-6, 
        2.584872967637478e-5, -2.6244994622929809e-4, 
        -0.00102234403086242668, 0.00863185023199066935, 
        0.01816610124346617164, -0.11151661591132254182, 
        -0.0897879180557114995, 0.23370422835726857838, 
        -2.957977086e-11, -6.150420279e-11, 
        4.88089472193e-9, 6.24218182859e-9, 
        -5.6561995267846e-7, -3.2077350100271e-7, 
        4.27856984896483e-5, -5.3069698280036e-7, 
        -0.0018714153878560023, 7.3172821883575326e-4, 
        0.03832854206202077374, -0.01874044416229564036, 
        -0.2302705940514488064, 0.05794254714300082167, 
        -1.203179184e-11, -3.5509887475e-10, 
        2.13620090322e-9, 4.813181550803e-8, 
        -2.6815627351126e-7, -4.39191902440609e-6, 
        2.211921535574542e-5, 2.45935054805247e-4, 
        -0.00106158930833741898, -0.00722666167492552731, 
        0.02396424370522244155, 0.08278313570069734843, 
        -0.15890724632166919294, -0.15383825653750118007
    };
    static double d[52] = {
        -1.272346002224188092e-14, 3.370464692346669075e-13, 
        -1.144940314335484869e-11, 6.863141561083429745e-10, 
        -9.491933932960924159e-8, 5.301676561445687562e-5, 
        0.162867503967639974, -3.652982212914147794e-13, 
        1.151126750560028914e-11, -5.165585095674343486e-10, 
        4.657991250060549892e-8, -1.186794704692706504e-5, 
        0.01562499999999994026, 
        -8.713069680903981555e-15, 3.140780373478474935e-13, 
        -1.139089186076256597e-11, 6.862299023338785566e-10, 
        -9.491926788274594674e-8, 5.301676558106268323e-5, 
        0.162867503967646622, -2.792555727162752006e-13, 
        1.108650207651756807e-11, -5.156745588549830981e-10, 
        4.657894859077370979e-8, -1.186794650130550256e-5, 
        0.01562499999987299901, 
        -6.304859171204770696e-15, 2.857249044208791652e-13, 
        -1.124956921556753188e-11, 6.858482894906716661e-10, 
        -9.49186795351689846e-8, 5.301676509057781574e-5, 
        0.1628675039678191167, -2.185193490132496053e-13, 
        1.048820673697426074e-11, -5.132819367467680132e-10, 
        4.65740943737299422e-8, -1.186794150862988921e-5, 
        0.01562499999779270706, 
        -4.74041720979200985e-15, 2.578715253644144182e-13, 
        -1.104148898414138857e-11, 6.850134201626289183e-10, 
        -9.49167823417491964e-8, 5.301676277588728159e-5, 
        0.1628675039690033136, -1.75512205749384229e-13, 
        9.848723331445182397e-12, -5.094535425482245697e-10, 
        4.656255982268609304e-8, -1.186792402114394891e-5, 
        0.01562499998712198636
    };

    if (x < 0.85) {
        t = x * x;
        y = (((((((a[0] * t + a[1]) * t + 
            a[2]) * t + a[3]) * t + a[4]) * t + 
            a[5]) * t + a[6]) * t + a[7]) * x;
        y = (((((((a[8] * t + a[9]) * t + 
            a[10]) * t + a[11]) * t + a[12]) * t + 
            a[13]) * t + a[14]) * t + a[15]) / x + 
            y * log(x);
    } else if (x < 4.15) {
        v = 4 / x;
        t = x - v;
        k = (int) (t + 4);
        t -= k - 3.5;
        k *= 14;
        y = (((((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12]) * t + b[k + 13]) * v;
    } else if (x < 12.5) {
        k = (int) x;
        t = x - (k + 0.5);
        k = 14 * (k - 4);
        y = ((((((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * t + c[k + 9]) * t + c[k + 10]) * t + 
            c[k + 11]) * t + c[k + 12]) * t + c[k + 13];
    } else {
        v = 24 / x;
        t = v * v;
        k = 13 * ((int) t);
        y = ((((((d[k] * t + d[k + 1]) * t + 
            d[k + 2]) * t + d[k + 3]) * t + d[k + 4]) * t + 
            d[k + 5]) * t + d[k + 6]) * sqrt(v);
        theta = (((((d[k + 7] * t + d[k + 8]) * t + 
            d[k + 9]) * t + d[k + 10]) * t + d[k + 11]) * t + 
            d[k + 12]) * v - 0.78539816339744830962;
        y *= -cos(x + theta);
    }
    return y;
}


double chebev(double a, double b, double c[], int m, double x)
{
  double d=0.0,dd=0.0,sv,y,y2;
  int j;
  
  if ((x-a)*(x-b) > 0.0) {
    fprintf(stderr,"x not in range in routine chebev\n");
    exit(1);
  }
  y2=2.0*(y=(2.0*x-a-b)/(b-a));
  for (j=m-1;j>=1;j--) {
    sv=d;
    d=y2*d-dd+c[j];
    dd=sv;
  }
  return y*d-dd+0.5*c[0];
}


#define NUSE1 5
#define NUSE2 5

void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi)
{
  double xx;
  static double c1[] = {
    -1.142022680371172e0,6.516511267076e-3,
    3.08709017308e-4,-3.470626964e-6,6.943764e-9,
    3.6780e-11,-1.36e-13};
  static double c2[] = {
    1.843740587300906e0,-0.076852840844786e0,
    1.271927136655e-3,-4.971736704e-6,-3.3126120e-8,
    2.42310e-10,-1.70e-13,-1.0e-15};
  
  xx=8.0*x*x-1.0;
  *gam1=chebev(-1.0,1.0,c1,NUSE1,xx);
  *gam2=chebev(-1.0,1.0,c2,NUSE2,xx);
  *gampl= *gam2-x*(*gam1);
  *gammi= *gam2+x*(*gam1);
}