bessy1.cpp

零阶、一阶和任意整数阶的第一二类贝赛尔函数算法

double bessy1(double x)
{
	double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
	double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,s7;
	double bbb,ccc,aaa,temp,xx,z,y;
    p1 = 1.0;                     p2 = 0.00183105;
    p3 = -0.00003516396496;       p4 = 0.000002457520174;
    p5 = -0.000000240337019;
    q1 = 0.04687499995;           q2 = -0.0002002690873;
    q3 = 0.000008449199096;       q4 = -0.00000088228987;
    q5 = 0.000000105787412;
    r1 = -4900604943000.0;        r2 = 1275274390000.0;
    r3 = -51534381390.0;          r4 = 734926455.1;
    r5 = -4237922.726;            r6 = 8511.937935;
    s1 = 24995805700000.0;        s2 = 424441966400.0;
    s3 = 3733650367.0;            s4 = 22459040.02;
    s5 = 102042.605;              s6 = 354.9632885;
    s7 = 1.0;
    if (x < 8.0)
	{
        y = x*x;
        aaa = r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))));
        bbb = s4 + y * (s5 + y * (s6 + y * s7));
        bbb = s1 + y * (s2 + y * (s3 + y * bbb));
        ccc = bessj1(x) * log(x) - 1.0 / x;
        temp = x * aaa / bbb + 0.636619772 * ccc;
	}
    else
	{
        z = 8.0 / x;
        y = z*z;
        xx = x - 2.356194491;
        aaa = sqrt(0.636619772 / x);
        bbb = p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)));
        ccc = q1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
        temp = aaa * (sin(xx) * bbb + z * cos(xx) * ccc);
	}
	return temp;
}