bessel.cpp

基于Visual C++环境的的Bessel函数程序

*************************************************************************/
double besseljvhankel(double n, double x)
{
    double result;
    double t;
    double u;
    double z;
    double k;
    double sg;
    double conv;
    double p;
    double q;
    double j;
    double m;
    double pp;
    double qq;
    int flag;

    m = 4.0*n*n;
    j = 1.0;
    z = 8.0*x;
    k = 1.0;
    p = 1.0;
    u = (m-1.0)/z;
    q = u;
    sg = 1.0;
    conv = 1.0;
    flag = 0;
    t = 1.0;
    pp = 1.0e38;
    qq = 1.0e38;
    while(t>ap::machineepsilon)
    {
        k = k+2.0;
        j = j+1.0;
        sg = -sg;
        u = u*((m-k*k)/(j*z));
        p = p+sg*u;
        k = k+2.0;
        j = j+1.0;
        u = u*((m-k*k)/(j*z));
        q = q+sg*u;
        t = fabs(u/p);
        if( t<conv )
        {
            conv = t;
            qq = q;
            pp = p;
            flag = 1;
        }
        if( flag!=0&&t>conv )
        {
            break;
        }
    }
    u = x-(0.5*n+0.25)*ap::pi();
    t = sqrt(2.0/(ap::pi()*x))*(pp*cos(u)-qq*sin(u));
    result = t;
    return result;
}


/*************************************************************************
Asymptotic expansion for large n.
AMS55 #9.3.35.
*************************************************************************/
double besseljnx(double n, double x)
{
    double result;
    double zeta;
    double sqz;
    double zz;
    double zp;
    double np;
    double cbn;
    double n23;
    double t;
    double z;
    double sz;
    double pp;
    double qq;
    double z32i;
    double zzi;
    double ak;
    double bk;
    double akl;
    double bkl;
    int sg;
    int doa;
    int dob;
    int nflg;
    int k;
    int s;
    int tk;
    int tkp1;
    int m;
    ap::real_1d_array u;
    double ai;
    double aip;
    double bi;
    double bip;
    ap::real_1d_array mu;
    ap::real_1d_array lambda;

    cbn = besseljvcbrt(n);
    z = (x-n)/cbn;
    if( fabs(z)<=0.7 )
    {
        result = besseljnt(n, x);
        return result;
    }
    z = x/n;
    zz = 1.0-z*z;
    if( zz==0.0 )
    {
        result = 0;
        return result;
    }
    if( zz>0.0 )
    {
        sz = sqrt(zz);
        t = 1.5*(log((1.0+sz)/z)-sz);
        zeta = pow(t*t, 1.0/3.0);
        nflg = 1;
    }
    else
    {
        sz = sqrt(-zz);
        t = 1.5*(sz-acos(1.0/z));
        zeta = -pow(t*t, 1.0/3.0);
        nflg = -1;
    }
    z32i = fabs(1.0/t);
    sqz = besseljvcbrt(t);
    n23 = besseljvcbrt(n*n);
    t = n23*zeta;
    airy(t, ai, aip, bi, bip);
    u.setbounds(0, 7);
    u(0) = 1.0;
    zzi = 1.0/zz;
    u(1) = -2.083333333333333333333333E-1;
    u(1) = u(1)*zzi+1.250000000000000000000000E-1;
    u(1) = u(1)/sz;
    u(2) = 3.342013888888888888888889E-1;
    u(2) = u(2)*zzi-4.010416666666666666666667E-1;
    u(2) = u(2)*zzi+7.031250000000000000000000E-2;
    u(2) = u(2)/zz;
    u(3) = -1.025812596450617283950617E+0;
    u(3) = u(3)*zzi+1.846462673611111111111111E+0;
    u(3) = u(3)*zzi-8.912109375000000000000000E-1;
    u(3) = u(3)*zzi+7.324218750000000000000000E-2;
    u(3) = u(3)/(sz*zz);
    pp = zz*zz;
    u(4) = 4.669584423426247427983539E+0;
    u(4) = u(4)*zzi-1.120700261622299382716049E+1;
    u(4) = u(4)*zzi+8.789123535156250000000000E+0;
    u(4) = u(4)*zzi-2.364086914062500000000000E+0;
    u(4) = u(4)*zzi+1.121520996093750000000000E-1;
    u(4) = u(4)/pp;
    u(5) = -2.8212072558200244877E1;
    u(5) = u(5)*zzi+8.4636217674600734632E1;
    u(5) = u(5)*zzi-9.1818241543240017361E1;
    u(5) = u(5)*zzi+4.2534998745388454861E1;
    u(5) = u(5)*zzi-7.3687943594796316964E0;
    u(5) = u(5)*zzi+2.27108001708984375E-1;
    u(5) = u(5)/(pp*sz);
    pp = pp*zz;
    u(6) = 2.1257013003921712286E2;
    u(6) = u(6)*zzi-7.6525246814118164230E2;
    u(6) = u(6)*zzi+1.0599904525279998779E3;
    u(6) = u(6)*zzi-6.9957962737613254123E2;
    u(6) = u(6)*zzi+2.1819051174421159048E2;
    u(6) = u(6)*zzi-2.6491430486951555525E1;
    u(6) = u(6)*zzi+5.7250142097473144531E-1;
    u(6) = u(6)/pp;
    u(7) = -1.9194576623184069963E3;
    u(7) = u(7)*zzi+8.0617221817373093845E3;
    u(7) = u(7)*zzi-1.3586550006434137439E4;
    u(7) = u(7)*zzi+1.1655393336864533248E4;
    u(7) = u(7)*zzi-5.3056469786134031084E3;
    u(7) = u(7)*zzi+1.2009029132163524628E3;
    u(7) = u(7)*zzi-1.0809091978839465550E2;
    u(7) = u(7)*zzi+1.7277275025844573975E0;
    u(7) = u(7)/(pp*sz);
    pp = 0.0;
    qq = 0.0;
    np = 1.0;
    mu.setbounds(0, 10);
    mu(0) = 1.0;
    mu(1) = -1.458333333333333333333333E-1;
    mu(2) = -9.874131944444444444444444E-2;
    mu(3) = -1.433120539158950617283951E-1;
    mu(4) = -3.172272026784135480967078E-1;
    mu(5) = -9.424291479571202491373028E-1;
    mu(6) = -3.511203040826354261542798E+0;
    mu(7) = -1.572726362036804512982712E+1;
    mu(8) = -8.228143909718594444224656E+1;
    mu(9) = -4.923553705236705240352022E+2;
    mu(10) = -3.316218568547972508762102E+3;
    lambda.setbounds(0, 10);
    lambda(0) = 1.0;
    lambda(1) = 1.041666666666666666666667E-1;
    lambda(2) = 8.355034722222222222222222E-2;
    lambda(3) = 1.282265745563271604938272E-1;
    lambda(4) = 2.918490264641404642489712E-1;
    lambda(5) = 8.816272674437576524187671E-1;
    lambda(6) = 3.321408281862767544702647E+0;
    lambda(7) = 1.499576298686255465867237E+1;
    lambda(8) = 7.892301301158651813848139E+1;
    lambda(9) = 4.744515388682643231611949E+2;
    lambda(10) = 3.207490090890661934704328E+3;
    doa = 1;
    dob = 1;
    akl = ap::maxrealnumber;
    bkl = ap::maxrealnumber;
    for(k = 0; k <= 3; k++)
    {
        tk = 2*k;
        tkp1 = tk+1;
        zp = 1.0;
        ak = 0.0;
        bk = 0.0;
        for(s = 0; s <= tk; s++)
        {
            if( doa!=0 )
            {
                if( s%4>1 )
                {
                    sg = nflg;
                }
                else
                {
                    sg = 1;
                }
                ak = ak+sg*mu(s)*zp*u(tk-s);
            }
            if( dob!=0 )
            {
                m = tkp1-s;
                if( (m+1)%4>1 )
                {
                    sg = nflg;
                }
                else
                {
                    sg = 1;
                }
                bk = bk+sg*lambda(s)*zp*u(m);
            }
            zp = zp*z32i;
        }
        if( doa!=0 )
        {
            ak = ak*np;
            t = fabs(ak);
            if( t<akl )
            {
                akl = t;
                pp = pp+ak;
            }
            else
            {
                doa = 0;
            }
        }
        if( dob!=0 )
        {
            bk = bk+lambda(tkp1)*zp*u(0);
            bk = -bk*np/sqz;
            t = fabs(bk);
            if( t<bkl )
            {
                bkl = t;
                qq = qq+bk;
            }
            else
            {
                dob = 0;
            }
        }
        if( np<ap::machineepsilon )
        {
            break;
        }
        np = np/n*n;
    }
    t = 4.0*zeta/zz;
    t = sqrt(sqrt(t));
    t = t*(ai*pp/besseljvcbrt(n)+aip*qq/(n23*n));
    result = t;
    return result;
}


/*************************************************************************
Asymptotic expansion for transition region,
n large and x close to n.
AMS55 #9.3.23.
*************************************************************************/
double besseljnt(double n, double x)
{
    double result;
    double z;
    double zz;
    double z3;
    double cbn;
    double n23;
    double cbtwo;
    double ai;
    double aip;
    double bi;
    double bip;
    double nk;
    double fk;
    double gk;
    double pp;
    double qq;
    int k;
    ap::real_1d_array f;
    ap::real_1d_array g;

    f.setbounds(0, 4);
    g.setbounds(0, 3);
    cbn = besseljvcbrt(n);
    z = (x-n)/cbn;
    cbtwo = besseljvcbrt(2.0);
    zz = -cbtwo*z;
    airy(zz, ai, aip, bi, bip);
    zz = z*z;
    z3 = zz*z;
    f(0) = 1.0;
    f(1) = -z/5.0;
    f(2) = -9.0000000000000000000e-2;
    f(2) = f(2)*z3+8.5714285714285714286e-2;
    f(2) = f(2)*zz;
    f(3) = 1.3671428571428571429e-1;
    f(3) = f(3)*z3-5.4920634920634920635e-2;
    f(3) = f(3)*z3-4.4444444444444444444e-3;
    f(4) = 1.3500000000000000000e-3;
    f(4) = f(4)*z3-1.6036054421768707483e-1;
    f(4) = f(4)*z3+4.2590187590187590188e-2;
    f(4) = f(4)*z3+2.7330447330447330447e-3;
    f(4) = f(4)*z;
    g(0) = 0.3*zz;
    g(1) = -2.4285714285714285714e-1;
    g(1) = g(1)*z3+1.4285714285714285714e-2;
    g(2) = -9.0000000000000000000e-3;
    g(2) = g(2)*z3+1.9396825396825396825e-1;
    g(2) = g(2)*z3-1.1746031746031746032e-2;
    g(2) = g(2)*z;
    g(3) = 1.9607142857142857143e-2;
    g(3) = g(3)*z3-1.5983694083694083694e-1;
    g(3) = g(3)*z3+6.3838383838383838384e-3;
    g(3) = g(3)*zz;
    pp = 0.0;
    qq = 0.0;
    nk = 1.0;
    n23 = besseljvcbrt(n*n);
    for(k = 0; k <= 4; k++)
    {
        fk = f(k)*nk;
        pp = pp+fk;
        if( k!=4 )
        {
            gk = g(k)*nk;
            qq = qq+gk;
        }
        nk = nk/n23;
    }
    fk = cbtwo*ai*pp/cbn+besseljvcbrt(4.0)*aip*qq/n;
    result = fk;
    return result;
}


/*************************************************************************
cube root
*************************************************************************/
double besseljvcbrt(double x)
{
    double result;

    if( x==0 )
    {
        result = 0;
    }
    if( x>0 )
    {
        result = pow(x, 1.0/3.0);
    }
    if( x<0 )
    {
        result = -pow(-x, 1.0/3.0);
    }
    return result;
}


/*************************************************************************
frexp
*************************************************************************/
double besselfrexp(double x, int& e)
{
    double result;
    double s;

    if( x==0 )
    {
        result = 0;
        e = 0;
        return result;
    }
    s = ap::sign(x);
    result = fabs(x);
    e = 0;
    while(result*1024*1024<1)
    {
        e = e-20;
        result = result*1024*1024;
    }
    while(result*1024<1)
    {
        e = e-10;
        result = result*1024;
    }
    while(result*2<1)
    {
        e = e-1;
        result = result*2;
    }
    while(result>=0.5*1024*1024)
    {
        e = e+20;
        result = result/(1024*1024);
    }
    while(result>=0.5*1024)
    {
        e = e+10;
        result = result/1024;
    }
    while(result>=0.5*2)
    {
        e = e