Return to j1.c CVS log

Up to [MW Coherent from dump] / coherent / b / lib / libm

Annotation of coherent/b/lib/libm/j1.c, revision 1.1

1.1     ! root        1: /*
        !             2:  * Evaluate the Bessel function of the first kind for the first order.
        !             3:  */
        !             4: 
        !             5: #include <math.h>
        !             6: 
        !             7: /*
        !             8:  * (Hart 6046, 17.13)
        !             9:  */
        !            10: static readonly double j1sntab[] ={
        !            11:         0.2214887880421963139207647803e+18,
        !            12:        -0.2512374214703212789513276482e+17,
        !            13:         0.8482420744781272654092270714e+15,
        !            14:        -0.1249820367262024853059386404e+14,
        !            15:         0.9316565296724673204941569090e+11,
        !            16:        -0.3686668987022981626057360734e+09,
        !            17:         0.7437023817119996441033971568e+06,
        !            18:        -0.6079530179607413599422162589e+03
        !            19: };
        !            20: static readonly double j1smtab[] ={
        !            21:         0.4429775760843926213068606431e+18,
        !            22:         0.5124712716484872112355190833e+16,
        !            23:         0.2989836307725487159974863805e+14,
        !            24:         0.1158192127466889329393351964e+12,
        !            25:         0.3281940344534196444499723799e+09,
        !            26:         0.6988586184485075744033771749e+06,
        !            27:         0.1077741289433304357312995618e+04,
        !            28:         0.1000000000000000000000000000e+01
        !            29: };
        !            30: 
        !            31: /*
        !            32:  * (Hart 6748, 16.96)
        !            33:  */
        !            34: static readonly double j1pntab[] ={
        !            35:        -0.1746957669440928588993119874e+07,
        !            36:        -0.2366943914052142834823993894e+07,
        !            37:        -0.8285036373872377453364829379e+06,
        !            38:        -0.7957456871350595920771321737e+05,
        !            39:        -0.1427906628827098109925825113e+04
        !            40: };
        !            41: static readonly double j1pmtab[] ={
        !            42:        -0.1746957669440928570099951216e+07,
        !            43:        -0.2363745139022653985015374700e+07,
        !            44:        -0.8242369906562818850190419776e+06,
        !            45:        -0.7814405008939111083484764129e+05,
        !            46:        -0.1308452983339079687508155034e+04,
        !            47:         0.1000000000000000000000000000e+01
        !            48: };
        !            49: 
        !            50: /*
        !            51:  * (Hart 7148, 17.16)
        !            52:  */
        !            53: static readonly double j1qntab[] ={
        !            54:         0.1437540512394930951685647580e+05,
        !            55:         0.2288127654437247527255486540e+05,
        !            56:         0.9758129726644802407557480600e+04,
        !            57:         0.1212097049668108451079609000e+04,
        !            58:         0.3134055238623227003160600000e+02
        !            59: };
        !            60: static readonly double j1qmtab[] ={
        !            61:         0.3066753093109186472663730271e+06,
        !            62:         0.4894441578927938446024906527e+06,
        !            63:         0.2102091447442188749228749142e+06,
        !            64:         0.2667395055960220714890160450e+05,
        !            65:         0.7531715785291537658631247000e+03,
        !            66:         0.1000000000000000000000000000e+01
        !            67: };
        !            68: 
        !            69: double
        !            70: j1(x)
        !            71: double x;
        !            72: {
        !            73:        double r, p, q, z, xn;
        !            74:        register int m;
        !            75: 
        !            76:        m = 0;
        !            77:        if (x < 0.0) {
        !            78:                x = -x;
        !            79:                m = 1;
        !            80:        }
        !            81:        if (x <= 8.0) {
        !            82:                r = x*x;
        !            83:                r = x*_pol(r, j1sntab, 8)/_pol(r, j1smtab, 8);
        !            84:        } else {
        !            85:                /* N.B. misprint in Hart 1968 edition, corrected 1978 reprint */
        !            86:                z = 8.0 / x;
        !            87:                xn = x - (3.0*PI)/4.0;
        !            88:                r = z*z;
        !            89:                p = _pol(r, j1pntab, 5) / _pol(r, j1pmtab, 6);
        !            90:                q = z*_pol(r, j1qntab, 5) / _pol(r, j1qmtab, 6);
        !            91:                r = sqrt(2.0/(PI*x)) * (p*cos(xn) - q*sin(xn));
        !            92:        }
        !            93:        if (m)
        !            94:                r = -r;
        !            95:        return (r);
        !            96: }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.