bessjy.cpp

Individual files are available in the following links: Bessjy -- Bessel functions Jn and Yn for r

    ecode = bessjy01a(x,jn[0],jn[1],yn[0],yn[1],jnp[0],jnp[1],ynp[0],ynp[1]);
    if (n < 2) return 0;
    bj0 = jn[0];
    bj1 = jn[1];
    if (n < (int)0.9*x) {
        for (k=2;k<=n;k++) {
            jn[k] = 2.0*(k-1.0)*bj1/x-bj0;
            bj0 = bj1;
            bj1 = jn[k];
        }
    }
    else {
        m = msta1(x,200);
        if (m < n) nm = m;
        else m = msta2(x,n,15);
        f2 = 0.0;
        f1 = 1.0e-100;
        for (k=m;k>=0;k--) {
            f = 2.0*(k+1.0)/x*f1-f2;
            if (k <= nm) jn[k] = f;
            f2 = f1;
            f1 = f;
        }
        if (fabs(bj0) > fabs(bj1)) cs = bj0/f;
        else cs = bj1/f2;
        for (k=0;k<=nm;k++) {
            jn[k] *= cs;
        }
    }
    for (k=2;k<=nm;k++) {
        jnp[k] = jn[k-1]-k*jn[k]/x;
    }
    f0 = yn[0];
    f1 = yn[1];
    for (k=2;k<=nm;k++) {
        f = 2.0*(k-1.0)*f1/x-f0;
        yn[k] = f;
        f0 = f1;
        f1 = f;
    }
    for (k=2;k<=nm;k++) {
        ynp[k] = yn[k-1]-k*yn[k]/x;
    }
    return 0;
}
//
//  Same input and output conventions as above. Different recurrence
//  relations used for 'x' < 300.
//
int bessjynb(int n,double x,int &nm,double *jn,double *yn,
    double *jnp,double *ynp)
{
    double t1,t2,f,f0,f1,f2,bj0,bj1,bjk,by0,by1,cu,s0,su,sv;
    double ec,bs,byk,p0,p1,q0,q1;
    static double a[] = {
        -0.7031250000000000e-1,
         0.1121520996093750,
        -0.5725014209747314,
         6.074042001273483};
    static double b[] = {
         0.7324218750000000e-1,
        -0.2271080017089844,
         1.727727502584457,
        -2.438052969955606e1};
    static double a1[] = {
         0.1171875,
        -0.1441955566406250,
         0.6765925884246826,
        -6.883914268109947};
    static double b1[] = {
       -0.1025390625,
        0.2775764465332031,
       -1.993531733751297,
        2.724882731126854e1};
        
    int i,k,m;
    nm = n;
    if ((x < 0.0) || (n < 0)) return 1;
    if (x < 1e-15) {
        for (i=0;i<=n;i++) {
            jn[i] = 0.0;
            yn[i] = -1e308;
            jnp[i] = 0.0;
            ynp[i] = 1e308;
        }
        jn[0] = 1.0;
        jnp[1] = 0.5;
        return 0;
    }
    if (x <= 300.0 || n > (int)(0.9*x)) {
        if (n == 0) nm = 1;
        m = msta1(x,200);
        if (m < nm) nm = m;
        else m = msta2(x,nm,15);
        bs = 0.0;
        su = 0.0;
        sv = 0.0;
        f2 = 0.0;
        f1 = 1.0e-100;
        for (k = m;k>=0;k--) {
            f = 2.0*(k+1.0)/x*f1 - f2;
            if (k <= nm) jn[k] = f;
            if ((k == 2*(int)(k/2)) && (k != 0)) {
                bs += 2.0*f;
//                su += pow(-1,k>>1)*f/(double)k;
                su += (-1)*((k & 2)-1)*f/(double)k;
            }
            else if (k > 1) {
//                sv += pow(-1,k>>1)*k*f/(k*k-1.0);
                sv += (-1)*((k & 2)-1)*(double)k*f/(k*k-1.0);
            }
            f2 = f1;
            f1 = f;
        }
        s0 = bs+f;
        for (k=0;k<=nm;k++) {
            jn[k] /= s0;
        }
        ec = log(0.5*x) +0.5772156649015329;
        by0 = M_2_PI*(ec*jn[0]-4.0*su/s0);
        yn[0] = by0;
        by1 = M_2_PI*((ec-1.0)*jn[1]-jn[0]/x-4.0*sv/s0);
        yn[1] = by1;
    }
    else {
        t1 = x-M_PI_4;
        p0 = 1.0;
        q0 = -0.125/x;
        for (k=0;k<4;k++) {
            p0 += a[k]*pow(x,-2*k-2);
            q0 += b[k]*pow(x,-2*k-3);
        }
        cu = sqrt(M_2_PI/x);
        bj0 = cu*(p0*cos(t1)-q0*sin(t1));
        by0 = cu*(p0*sin(t1)+q0*cos(t1));
        jn[0] = bj0;
        yn[0] = by0;
        t2 = x-0.75*M_PI;
        p1 = 1.0;
        q1 = 0.375/x;
        for (k=0;k<4;k++) {
            p1 += a1[k]*pow(x,-2*k-2);
            q1 += b1[k]*pow(x,-2*k-3);
        }
        bj1 = cu*(p1*cos(t2)-q1*sin(t2));
        by1 = cu*(p1*sin(t2)+q1*cos(t2));
        jn[1] = bj1;
        yn[1] = by1;
        for (k=2;k<=nm;k++) {
            bjk = 2.0*(k-1.0)*bj1/x-bj0;
            jn[k] = bjk;
            bj0 = bj1;
            bj1 = bjk;
        }
    }
    jnp[0] = -jn[1];
    for (k=1;k<=nm;k++) {
        jnp[k] = jn[k-1]-k*jn[k]/x;
    }
    for (k=2;k<=nm;k++) {
        byk = 2.0*(k-1.0)*by1/x-by0;
        yn[k] = byk;
        by0 = by1;
        by1 = byk;
    }
    ynp[0] = -yn[1];
    for (k=1;k<=nm;k++) {
        ynp[k] = yn[k-1]-k*yn[k]/x;
    }
    return 0;

}

//  The following routine computes Bessel Jv(x) and Yv(x) for
//  arbitrary positive order (v). For negative order, use:
//
//      J-v(x) = Jv(x)cos(v pi) - Yv(x)sin(v pi)
//      Y-v(x) = Jv(x)sin(v pi) + Yv(x)cos(v pi)
//
int bessjyv(double v,double x,double &vm,double *jv,double *yv,
    double *djv,double *dyv)
{
    double v0,vl,vg,vv,a,a0,r,x2,bjv0,bjv1,bjvl,f,f0,f1,f2;
    double r0,r1,ck,cs,cs0,cs1,sk,qx,px,byv0,byv1,rp,xk,rq;
    double b,ec,w0,w1,bjy0,bjy1,bju0,bju1,pv0,pv1,byvk;
    int j,k,l,m,n,kz;

    x2 = x*x;
    n = (int)v;
    v0 = v-n;
    if ((x < 0.0) || (v < 0.0)) return 1;
    if (x < 1e-15) {
        for (k=0;k<=n;k++) {
            jv[k] = 0.0;
            yv[k] = -1e308;
            djv[k] = 0.0;
            dyv[k] = 1e308;
            if (v0 == 0.0) {
                jv[0] = 1.0;
                djv[1] = 0.5;
            }
            else djv[0] = 1e308;
        }
        vm = v;
        return 0;
    }
    if (x <= 12.0) {
        for (l=0;l<2;l++) {
            vl = v0 + l;
            bjvl = 1.0;
            r = 1.0;
            for (k=1;k<=40;k++) {
                r *= -0.25*x2/(k*(k+vl));
                bjvl += r;
                if (fabs(r) < fabs(bjvl)*1e-15) break;
            }
            vg = 1.0 + vl;
            a = pow(0.5*x,vl)/gamma(vg);
            if (l == 0) bjv0 = bjvl*a;
            else bjv1 = bjvl*a;
        }
    }
    else {
        if (x >= 50.0) kz = 8;
        else if (x >= 35.0) kz = 10;
        else kz = 11;
        for (j=0;j<2;j++) {
            vv = 4.0*(j+v0)*(j+v0);
            px = 1.0;
            rp = 1.0;
            for (k=1;k<=kz;k++) {
                rp *= (-0.78125e-2)*(vv-pow(4.0*k-3.0,2.0))*
                    (vv-pow(4.0*k-1.0,2.0))/(k*(2.0*k-1.0)*x2);
                px += rp;
            }
            qx = 1.0;
            rq = 1.0;
            for (k=1;k<=kz;k++) {
                rq *= (-0.78125e-2)*(vv-pow(4.0*k-1.0,2.0))*
                    (vv-pow(4.0*k+1.0,2.0))/(k*(2.0*k+1.0)*x2);
                qx += rq;
            }
            qx *= 0.125*(vv-1.0)/x;
            xk = x-(0.5*(j+v0)+0.25)*M_PI;
            a0 = sqrt(M_2_PI/x);
            ck = cos(xk);
            sk = sin(xk);

            if (j == 0) {
                bjv0 = a0*(px*ck-qx*sk);
                byv0 = a0*(px*sk+qx*ck);
            }
            else {
                bjv1 = a0*(px*ck-qx*sk);
                byv1 = a0*(px*sk+qx*ck);
            }
        }
    }
    jv[0] = bjv0;
    jv[1] = bjv1;
    djv[0] = v0*jv[0]/x-jv[1];
    djv[1] = -(1.0+v0)*jv[1]/x+jv[0];
    if ((n >= 2) && (n <= (int)(0.9*x))) {
        f0 = bjv0;
        f1 = bjv1;
        for (k=2;k<=n;k++) {
            f = 2.0*(k+v0-1.0)*f1/x-f0;
            jv[k] = f;
            f0 = f1;
            f1 = f;
        }
    }
    else if (n >= 2) {
        m = msta1(x,200);
        if (m < n) n = m;
        else m = msta2(x,n,15);
        f2 = 0.0;
        f1 = 1.0e-100;
        for (k=m;k>=0;k--) {
            f = 2.0*(v0+k+1.0)/x*f1-f2;
            if (k <= n) jv[k] = f;
            f2 = f1;
            f1 = f;
        }
        if (fabs(bjv0) > fabs(bjv1)) cs = bjv0/f;
        else cs = bjv1/f2;
        for (k=0;k<=n;k++) {
            jv[k] *= cs;
        }
    }
    for (k=2;k<=n;k++) {
        djv[k] = -(k+v0)*jv[k]/x+jv[k-1];
    }
    if (x <= 12.0) {
        if (v0 != 0.0) {
            for (l=0;l<2;l++) {
                vl = v0 +l;
                bjvl = 1.0;
                r = 1.0;
                for (k=1;k<=40;k++) {
                    r *= -0.25*x2/(k*(k-vl));
                    bjvl += r;
                    if (fabs(r) < fabs(bjvl)*1e-15) break;
                }
                vg = 1.0-vl;
                b = pow(2.0/x,vl)/gamma(vg);
                if (l == 0) bju0 = bjvl*b;
                else bju1 = bjvl*b;
            }
            pv0 = M_PI*v0;
            pv1 = M_PI*(1.0+v0);
            byv0 = (bjv0*cos(pv0)-bju0)/sin(pv0);
            byv1 = (bjv1*cos(pv1)-bju1)/sin(pv1);
        }
        else {
            ec = log(0.5*x)+el;
            cs0 = 0.0;
            w0 = 0.0;
            r0 = 1.0;
            for (k=1;k<=30;k++) {
                w0 += 1.0/k;
                r0 *= -0.25*x2/(k*k);
                cs0 += r0*w0;
            }
            byv0 = M_2_PI*(ec*bjv0-cs0);
            cs1 = 1.0;
            w1 = 0.0;
            r1 = 1.0;
            for (k=1;k<=30;k++) {
                w1 += 1.0/k;
                r1 *= -0.25*x2/(k*(k+1));
                cs1 += r1*(2.0*w1+1.0/(k+1.0));
            }
            byv1 = M_2_PI*(ec*bjv1-1.0/x-0.25*x*cs1);
        }
    }
    yv[0] = byv0;
    yv[1] = byv1;
    for (k=2;k<=n;k++) {
        byvk = 2.0*(v0+k-1.0)*byv1/x-byv0;
        yv[k] = byvk;
        byv0 = byv1;
        byv1 = byvk;
    }
    dyv[0] = v0*yv[0]/x-yv[1];
    for (k=1;k<=n;k++) {
        dyv[k] = -(k+v0)*yv[k]/x+yv[k-1];
    }
    vm = n + v0;
    return 0;
}