cbessjy.cpp

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

    }
    if ((a0 <= 300.0) || (n > (int)(0.25*a0))) {
        if (n == 0) nm = 1;
        m = msta1(a0,200);
        if (m < nm) nm = m;
        else m = msta2(a0,nm,15);
        cbs = czero;
        csu = czero;
        csv = czero;
        cf2 = czero;
        cf1 = complex<double> (1.0e-100,0.0);
        for (k=m;k>=0;k--) {
            cf = 2.0*(k+1.0)*cf1/z-cf2;
            if (k <= nm) cj[k] = cf;
            if (((k & 1) == 0) && (k != 0)) {
                if (y0 <= 1.0) {
                    cbs += 2.0*cf;
                }
                else {
                    cbs += (-1)*((k & 2)-1)*2.0*cf;
                }
                csu += (double)((-1)*((k & 2)-1))*cf/(double)k;
            }
            else if (k > 1) {
                csv += (double)((-1)*((k & 2)-1)*k)*cf/(double)(k*k-1.0);    
            }
            cf2 = cf1;
            cf1 = cf;
        }
        if (y0 <= 1.0) cs0 = cbs+cf;
        else cs0 = (cbs+cf)/cos(z);
        for (k=0;k<=nm;k++) {
            cj[k] /= cs0;
        }
        ce = log(0.5*z)+el;
        cy[0] = M_2_PI*(ce*cj[0]-4.0*csu/cs0);
        cy[1] = M_2_PI*(-cj[0]/z+(ce-1.0)*cj[1]-4.0*csv/cs0);
    }
    else {
        ct1 = z-M_PI_4;
        cp0 = cone;
        for (k=0;k<4;k++) {
            cp0 += a[k]*pow(z,-2.0*k-2.0);
        }
        cq0 = -0.125/z;
        for (k=0;k<4;k++) {
            cq0 += b[k] *pow(z,-2.0*k-3.0);
        }
        cu = sqrt(M_2_PI/z);
        cbj0 = cu*(cp0*cos(ct1)-cq0*sin(ct1));
        cby0 = cu*(cp0*sin(ct1)+cq0*cos(ct1));
        cj[0] = cbj0;
        cy[0] = cby0;
        ct2 = z-0.75*M_PI;
        cp1 = cone;
        for (k=0;k<4;k++) {
            cp1 += a1[k]*pow(z,-2.0*k-2.0);
        }
        cq1 = 0.375/z;
        for (k=0;k<4;k++) {
            cq1 += b1[k]*pow(z,-2.0*k-3.0);
        }
        cbj1 = cu*(cp1*cos(ct2)-cq1*sin(ct2));
        cby1 = cu*(cp1*sin(ct2)+cq1*sin(ct2));
        cj[1] = cbj1;
        cy[1] = cby1;
        for (k=2;k<=n;k++) {
            cbjk = 2.0*(k-1.0)*cbj1/z-cbj0;
            cj[k] = cbjk;
            cbj0 = cbj1;
            cbj1 = cbjk;
        }
    }
    cjp[0] = -cj[1];
    for (k=1;k<=nm;k++) {
        cjp[k] = cj[k-1]-(double)k*cj[k]/z;
    }
    if (abs(cj[0]) > 1.0)
        cy[1] = (cj[1]*cy[0]-2.0/(M_PI*z))/cj[0];
    for (k=2;k<=nm;k++) {
        if (abs(cj[k-1]) >= abs(cj[k-2]))
            cyy = (cj[k]*cy[k-1]-2.0/(M_PI*z))/cj[k-1];
        else
            cyy = (cj[k]*cy[k-2]-4.0*(k-1.0)/(M_PI*z*z))/cj[k-2];
        cy[k] = cyy;
    }
    cyp[0] = -cy[1];
    for (k=1;k<=nm;k++) {
        cyp[k] = cy[k-1]-(double)k*cy[k]/z;
    }

    return 0;
}

int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv,
    complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp)
{
    complex<double> z1,z2,zk,cjvl,cr,ca,cjv0,cjv1,cpz,crp;
    complex<double> cqz,crq,ca0,cck,csk,cyv0,cyv1,cju0,cju1,cb;
    complex<double> cs,cs0,cr0,cs1,cr1,cec,cf,cf0,cf1,cf2;
    complex<double> cfac0,cfac1,cg0,cg1,cyk,cp11,cp12,cp21,cp22;
    complex<double> ch0,ch1,ch2,cyl1,cyl2,cylk;

    double a0,v0,pv0,pv1,vl,ga,gb,vg,vv,w0,w1,ya0,yak,ya1,wa;
    int j,n,k,kz,l,lb,lb0,m;

    a0 = abs(z);
    z1 = z;
    z2 = z*z;
    n = (int)v;


    v0 = v-n;

    pv0 = M_PI*v0;
    pv1 = M_PI*(1.0+v0);
    if (a0 < 1.0e-100) {
        for (k=0;k<=n;k++) {
            cjv[k] = czero;
            cyv[k] = complex<double> (-1e308,0);
            cjvp[k] = czero;
            cyvp[k] = complex<double> (1e308,0);

        }
        if (v0 == 0.0) {
            cjv[0] = cone;
            cjvp[1] = complex<double> (0.5,0.0);
        }
        else {
            cjvp[0] = complex<double> (1e308,0);
        }
        vm = v;
        return 0;
    }
    if (real(z1) < 0.0) z1 = -z;
    if (a0 <= 12.0) {
        for (l=0;l<2;l++) {
            vl = v0+l;
            cjvl = cone;
            cr = cone;
            for (k=1;k<=40;k++) {
                cr *= -0.25*z2/(k*(k+vl));
                cjvl += cr;
                if (abs(cr) < abs(cjvl)*eps) break;
            }
           vg = 1.0 + vl;
           ga = gamma(vg);
           ca = pow(0.5*z1,vl)/ga;
           if (l == 0) cjv0 = cjvl*ca;
           else cjv1 = cjvl*ca;
        }
    }
    else {
        if (a0 >= 50.0) kz = 8;
        else if (a0 >= 35.0) kz = 10;
        else kz = 11;
        for (j=0;j<2;j++) {
            vv = 4.0*(j+v0)*(j+v0);
            cpz = cone;
            crp = cone;
            for (k=1;k<=kz;k++) {
                crp = -0.78125e-2*crp*(vv-pow(4.0*k-3.0,2.0))*
                    (vv-pow(4.0*k-1.0,2.0))/(k*(2.0*k-1.0)*z2);
                cpz += crp;
            }
            cqz = cone;
            crq = cone;
            for (k=1;k<=kz;k++) {
                crq = -0.78125e-2*crq*(vv-pow(4.0*k-1.0,2.0))*
                    (vv-pow(4.0*k+1.0,2.0))/(k*(2.0*k+1.0)*z2);
                cqz += crq;
            }
            cqz *= 0.125*(vv-1.0)/z1;
            zk = z1-(0.5*(j+v0)+0.25)*M_PI;
            ca0 = sqrt(M_2_PI/z1);
            cck = cos(zk);
            csk = sin(zk);
            if (j == 0) {
                cjv0 = ca0*(cpz*cck-cqz*csk);
                cyv0 = ca0*(cpz*csk+cqz+cck);
            }
            else {
                cjv1 = ca0*(cpz*cck-cqz*csk);
                cyv1 = ca0*(cpz*csk+cqz*cck);
            }
        }
    }
    if (a0 <= 12.0) {
        if (v0 != 0.0) {
            for (l=0;l<2;l++) {
                vl = v0+l;
                cjvl = cone;
                cr = cone;
                for (k=1;k<=40;k++) {
                    cr *= -0.25*z2/(k*(k-vl));
                    cjvl += cr;
                    if (abs(cr) < abs(cjvl)*eps) break;
                }
                vg = 1.0-vl;
                gb = gamma(vg);
                cb = pow(2.0/z1,vl)/gb;
                if (l == 0) cju0 = cjvl*cb;
                else cju1 = cjvl*cb;
            }
            cyv0 = (cjv0*cos(pv0)-cju0)/sin(pv0);
            cyv1 = (cjv1*cos(pv1)-cju1)/sin(pv1);
        }
        else {
            cec = log(0.5*z1)+el;
            cs0 = czero;
            w0 = 0.0;
            cr0 = cone;
            for (k=1;k<=30;k++) {
                w0 += 1.0/k;
                cr0 *= -0.25*z2/(double)(k*k);
                cs0 += cr0*w0;
            }
            cyv0 = M_2_PI*(cec*cjv0-cs0);
            cs1 = cone;
            w1 = 0.0;
            cr1 = cone;
            for (k=1;k<=30;k++) {
                w1 += 1.0/k;
                cr1 *= -0.25*z2/(k*(k+1.0));
                cs1 += cr1*(2.0*w1+1.0/(k+1.0));
            }
            cyv1 = M_2_PI*(cec*cjv1-1.0/z1-0.25*z1*cs1);
        }
    }
    if (real(z) < 0.0) {
        cfac0 = exp(pv0*cii);
        cfac1 = exp(pv1*cii);
        if (imag(z) < 0.0) {
            cyv0 = cfac0*cyv0-2.0*cii*cos(pv0)*cjv0;
            cyv1 = cfac1*cyv1-2.0*cii*cos(pv1)*cjv1;
            cjv0 /= cfac0;
            cjv1 /= cfac1;
        }
        else if (imag(z) > 0.0) {
            cyv0 = cyv0/cfac0+2.0*cii*cos(pv0)*cjv0;
            cyv1 = cyv1/cfac1+2.0*cii*cos(pv1)*cjv1;
            cjv0 *= cfac0;
            cjv1 *= cfac1;
        }
    }
    cjv[0] = cjv0;
    cjv[1] = cjv1;
    if ((n >= 2) && (n <= (int)(0.25*a0))) {
        cf0 = cjv0;
        cf1 = cjv1;
        for (k=2;k<= n;k++) {
            cf = 2.0*(k+v0-1.0)*cf1/z-cf0;
            cjv[k] = cf;
            cf0 = cf1;
            cf1 = cf;
        }
    }
    else if (n >= 2) {
        m = msta1(a0,200);
        if (m < n) n = m;
        else  m = msta2(a0,n,15);
        cf2 = czero;
        cf1 = complex<double>(1.0e-100,0.0);
        for (k=m;k>=0;k--) {
            cf = 2.0*(v0+k+1.0)*cf1/z-cf2;
            if (k <= n) cjv[k] = cf;
            cf2 = cf1;
            cf1 = cf;
        }
        if (abs(cjv0) > abs(cjv1)) cs = cjv0/cf;
        else cs = cjv1/cf2;
        for (k=0;k<=n;k++) {
            cjv[k] *= cs;
        }
    }
    cjvp[0] = v0*cjv[0]/z-cjv[1];
    for (k=1;k<=n;k++) {
        cjvp[k] = -(k+v0)*cjv[k]/z+cjv[k-1];
    }
    cyv[0] = cyv0;
    cyv[1] = cyv1;
    ya0 = abs(cyv0);
    lb = 0;
    cg0 = cyv0;
    cg1 = cyv1;
    for (k=2;k<=n;k++) {
        cyk = 2.0*(v0+k-1.0)*cg1/z-cg0;
        yak = abs(cyk);
        ya1 = abs(cg0);
        if ((yak < ya0) && (yak< ya1)) lb = k;
        cyv[k] = cyk;
        cg0 = cg1;
        cg1 = cyk;
    }
    lb0 = 0;
    if ((lb > 4) && (imag(z) != 0.0)) {
        while(lb != lb0) {
            ch2 = cone;
            ch1 = czero;
            lb0 = lb;
            for (k=lb;k>=1;k--) {
                ch0 = 2.0*(k+v0)*ch1/z-ch2;
                ch2 = ch1;
                ch1 = ch0;
            }
            cp12 = ch0;
            cp22 = ch2;
            ch2 = czero;
            ch1 = cone;
            for (k=lb;k>=1;k--) {
                ch0 = 2.0*(k+v0)*ch1/z-ch2;
                ch2 = ch1;
                ch1 = ch0;
            }
            cp11 = ch0;
            cp21 = ch2;
            if (lb == n)
                cjv[lb+1] = 2.0*(lb+v0)*cjv[lb]/z-cjv[lb-1];
            if (abs(cjv[0]) > abs(cjv[1])) {
                cyv[lb+1] = (cjv[lb+1]*cyv0-2.0*cp11/(M_PI*z))/cjv[0];
                cyv[lb] = (cjv[lb]*cyv0+2.0*cp12/(M_PI*z))/cjv[0];
            }
            else {
                cyv[lb+1] = (cjv[lb+1]*cyv1-2.0*cp21/(M_PI*z))/cjv[1];
                cyv[lb] = (cjv[lb]*cyv1+2.0*cp22/(M_PI*z))/cjv[1];
            }
            cyl2 = cyv[lb+1];
            cyl1 = cyv[lb];
            for (k=lb-1;k>=0;k--) {
                cylk = 2.0*(k+v0+1.0)*cyl1/z-cyl2;
                cyv[k] = cylk;
                cyl2 = cyl1;
                cyl1 = cylk;
            }
            cyl1 = cyv[lb];
            cyl2 = cyv[lb+1];
            for (k=lb+1;k<n;k++) {
                cylk = 2.0*(k+v0)*cyl2/z-cyl1;
                cyv[k+1] = cylk;
                cyl1 = cyl2;
                cyl2 = cylk;
            }
            for (k=2;k<=n;k++) {
                wa = abs(cyv[k]);
                if (wa < abs(cyv[k-1])) lb = k;
            }
        }
    }
    cyvp[0] = v0*cyv[0]/z-cyv[1];
    for (k=1;k<=n;k++) {
        cyvp[k] = cyv[k-1]-(k+v0)*cyv[k]/z;
    }
    vm = n+v0;
    return 0;
}