double bessk1(double x)
{
	double p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7;
	double y,aaa,bbb,ccc,temp; 
    p1 = 1.0;               p2 = 0.15443144;
    p3 = -0.67278579;       p4 = -0.18156897;
    p5 = -0.01919402;       p6 = -0.00110404;
    p7 = -0.00004686;
    q1 = 1.25331414;        q2 = 0.23498619;
    q3 = -0.0365562;        q4 = 0.01504268;
    q5 = -0.00780353;       q6 = 0.00325614;
    q7 = -0.00068245;
    if (x <= 2.0)
	{
        y = x * x / 4;
        aaa = log(x / 2.0) * bessi1(x);
        ccc = y * (p5 + y * (p6 + y * p7));
        bbb = p1 + y * (p2 + y * (p3 + y * (p4 + ccc)));
        temp = aaa + (1.0 / x) * bbb;
	}
    else
	{
        y = (2.0 / x);
        bbb = y * (q5 + y * (q6 + y * q7));
        aaa = q1 + y * (q2 + y * (q3 + y * (q4 + bbb)));
        temp = (exp(-x) / sqrt(x)) * aaa;
	}    
	return temp;
}