double bessi1(double x)
{
	double p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9;
	double y,aaa,bbb,ax,temp;
    p1 = 0.5;              p2 = 0.87890594;
    p3 = 0.51498869;       p4 = 0.15084934;
    p5 = 0.02658733;       p6 = 0.00301532;
    p7 = 0.00032411;
    q1 = 0.39894228;       q2 = -0.03988024;
    q3 = -0.00362018;      q4 = 0.00163801;
    q5 = -0.01031555;      q6 = 0.02282967;
    q7 = -0.02895312;      q8 = 0.01787654;
    q9 = -0.00420059;
    if (fabs(x) < 3.75 )
	{
        y = (x / 3.75) *(x / 3.75) ;
        aaa = y * (p4 + y * (p5 + y * (p6 + y * p7)));
        temp = x * (p1 + y * (p2 + y * (p3 + aaa)));
	}
    else
	{
        ax = fabs(x);
        y = 3.75 / ax;
        aaa = exp(ax) / sqrt(ax);
        bbb = y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))));
        temp = aaa * (q1 + y * (q2 + y * (q3 + y * (q4 + bbb))));
	}
    return temp;
}