Function BESSI(N:integer; X:real):real;
const
    IACC = 40;    BIGNO = 1.e10;
    BIGNI = 1.e-10;
var
    TOX,BIP,BI,BIM,ANS:real; M,J:integer;
begin
    If N < 2 Then
        ShowMessage('bad argument N in BESSI');
    TOX:=2 / X;
    BIP:=0;
    BI:=1;
    BESSI:=0;
    M:=2 * ((N + Trunc(Sqrt(IACC * N))));
    For J:=M DownTo 1 do
    begin
        BIM:=BIP + J * TOX * BI;
        BIP:=BI;
        BI:=BIM;
        If Abs(BI) > BIGNO Then
        begin
            ANS:=ANS * BIGNI;
            BI:=BI * BIGNI;
            BIP:=BIP * BIGNI;
        end;
        If J = N Then ANS:=BIP;
    end;
    BESSI:=ANS * BESSI0(X) / BI;
end;