unit Unit2;

interface
Function BESSI1(X:Real):real;

implementation
Function BESSI1(X:Real):real;
const
    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;
var
    Y,AAA,AX,BBB:Real;
begin
    If Abs(X) < 3.75 Then
    begin
        Y:=Sqr(X / 3.75);
        AAA:=Y * (P4 + Y * (P5 + Y * (P6 + Y * P7)));
        BESSI1:=X * (P1 + Y * (P2 + Y * (P3 + AAA)));
    end
    Else
    begin
        AX:=Abs(X);
        Y:=3.75 / AX;
        AAA:=Exp(AX) / Sqrt(AX);
        BBB:=Y * (Q5 + Y * (Q6 + Y * (Q7 + Y * (Q8 + Y * Q9))));
        BESSI1:=AAA * (Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + BBB))));
    end;
end;

end.