unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs;
Function BESSI0(X:REAL):REAL;
Function BESSI(N:integer; X:real):real;

implementation
Function BESSI0(X:REAL):REAL;
CONST
    P1 = 1;             P2 = 3.5156229;
    P3 = 3.0899424;     P4 = 1.2067492;
    P5 = 0.2659732;     P6 = 0.0360768;
    P7 = 0.0045813;
    Q1 = 0.39894228;    Q2 = 0.01328592;
    Q3 = 0.00225319;    Q4 = -0.00157565;
    Q5 = 0.00916281;    Q6 = -0.02057706;
    Q7 = 0.02635537;    Q8 = -0.01647633;
    Q9 = 0.00392377;
VAR
    Y,AAA,BBB,AX:Real;
begin
    If Abs(X) < 3.75 Then
      begin
        Y:=Sqr(X / 3.75);
        AAA:=Y * (P5 + Y * (P6 + Y * P7));
        BESSI0:=P1 + Y * (P2 + Y * (P3 + Y * (P4 + AAA)));
      end
    Else
      begin
        AX:=Abs(X);
        Y:=3.75 / AX;
        AAA:=Exp(AX) / Sqrt(AX);
        BBB:=Q4 + Y * (Q5 + Y * (Q6 + Y * (Q7 + Y * (Q8 + Y * Q9))));
        BESSI0:=AAA * (Q1 + Y * (Q2 + Y * (Q3 + Y * BBB)));
      end;
end;

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;

end.