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

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.