Return to besjn.c CVS log

Up to [Research Unix] / researchv10dc / libc / math

Annotation of researchv10dc/libc/math/besjn.c, revision 1.1

1.1     ! root        1: /*
        !             2:        floating point Bessel's function of
        !             3:        the first and second kinds and of
        !             4:        integer order.
        !             5: 
        !             6:        int n;
        !             7:        double x;
        !             8:        besjn(n,x);
        !             9: 
        !            10:        returns the value of Jn(x) for all
        !            11:        integer values of n and all real values
        !            12:        of x.
        !            13: 
        !            14:        There are no error returns.
        !            15:        Calls besj0, besj1.
        !            16: 
        !            17:        For n=0, besj0(x) is called,
        !            18:        for n=1, besj1(x) is called,
        !            19:        for n<x, forward recursion us used starting
        !            20:        from values of besj0(x) and besj1(x).
        !            21:        for n>x, a continued fraction approximation to
        !            22:        j(n,x)/j(n-1,x) is evaluated and then backward
        !            23:        recursion is used starting from a supposed value
        !            24:        for j(n,x). The resulting value of j(0,x) is
        !            25:        compared with the actual value to correct the
        !            26:        supposed value of j(n,x).
        !            27: 
        !            28:        yn(n,x) is similar in all respects, except
        !            29:        that forward recursion is used for all
        !            30:        values of n>1.
        !            31: */
        !            32: 
        !            33: #include <math.h>
        !            34: #include <errno.h>
        !            35: 
        !            36: int    errno;
        !            37: 
        !            38: double
        !            39: besjn(n,x) int n; double x;{
        !            40:        int i;
        !            41:        double a, b, temp;
        !            42:        double xsq, t;
        !            43:        double besj0(), besj1();
        !            44: 
        !            45:        if(n<0){
        !            46:                n = -n;
        !            47:                x = -x;
        !            48:        }
        !            49:        if(n==0) return(besj0(x));
        !            50:        if(n==1) return(besj1(x));
        !            51:        if(x == 0.) return(0.);
        !            52:        if(n>x) goto recurs;
        !            53: 
        !            54:        a = besj0(x);
        !            55:        b = besj1(x);
        !            56:        for(i=1;i<n;i++){
        !            57:                temp = b;
        !            58:                b = (2.*i/x)*b - a;
        !            59:                a = temp;
        !            60:        }
        !            61:        return(b);
        !            62: 
        !            63: recurs:
        !            64:        xsq = x*x;
        !            65:        for(t=0,i=n+16;i>n;i--){
        !            66:                t = xsq/(2.*i - t);
        !            67:        }
        !            68:        t = x/(2.*n-t);
        !            69: 
        !            70:        a = t;
        !            71:        b = 1;
        !            72:        for(i=n-1;i>0;i--){
        !            73:                temp = b;
        !            74:                b = (2.*i/x)*b - a;
        !            75:                a = temp;
        !            76:        }
        !            77:        return(t*besj0(x)/b);
        !            78: }
        !            79: 
        !            80: double
        !            81: besyn(n,x) int n; double x;{
        !            82:        int i;
        !            83:        int sign;
        !            84:        double a, b, temp;
        !            85:        double besy0(), besy1();
        !            86: 
        !            87:        if (x <= 0) {
        !            88:                errno = EDOM;
        !            89:                return(-HUGE);
        !            90:        }
        !            91:        sign = 1;
        !            92:        if(n<0){
        !            93:                n = -n;
        !            94:                if(n%2 == 1) sign = -1;
        !            95:        }
        !            96:        if(n==0) return(besy0(x));
        !            97:        if(n==1) return(sign*besy1(x));
        !            98: 
        !            99:        a = besy0(x);
        !           100:        b = besy1(x);
        !           101:        for(i=1;i<n;i++){
        !           102:                temp = b;
        !           103:                b = (2.*i/x)*b - a;
        !           104:                a = temp;
        !           105:        }
        !           106:        return(sign*b);
        !           107: }

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