unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs;
Function BESSI0(X:REAL):REAL;
Function BESSI1(X:Real):real;
Function BESSK1(X:real):real;
Function BESSK0(X:Real):Real;
Function BESSK(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 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;

Function BESSK1(X:real):real;
const
    P1 = 1;    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;
var
    AAA,BBB,CCC,Y:REAL;
BEGIN
    If X <= 2 Then
      BEGIN
        Y:=X * X / 4;
        AAA:=Ln(X / 2) * BESSI1(X);
        CCC:=Y * (P5 + Y * (P6 + Y * P7));
        BBB:=P1 + Y * (P2 + Y * (P3 + Y * (P4 + CCC)));
        BESSK1:=AAA + (1 / X) * BBB;
      end
    Else
      begin
        Y:=(2 / X);
        BBB:=Y * (Q5 + Y * (Q6 + Y * Q7));
        AAA:=Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + BBB)));
        BESSK1:=(Exp(-X) / Sqrt(X)) * AAA;
      end;
end;

Function BESSK0(X:Real):Real;
const
    P1 = -0.57721566;    P2 = 0.4227842;
    P3 = 0.23069756;    P4 = 0.0348859;
    P5 = 0.00262698;    P6 = 0.0001075;
    P7 = 0.0000074;
    Q1 = 1.25331414;    Q2 = -0.07832358;
    Q3 = 0.02189568;    Q4 = -0.01062446;
    Q5 = 0.00587872;    Q6 = -0.0025154;
    Q7 = 0.00053208;
var
    Y,AAA,BBB:Real;
begin
    If X <= 2 Then
      begin
        Y:=X * X / 4;
        BBB:=Y * (P5 + Y * (P6 + Y * P7));
        AAA:=P1 + Y * (P2 + Y * (P3 + Y * (P4 + BBB)));
        BESSK0:=(-Ln(X / 2) * BESSI0(X)) + AAA;
      end
    Else
      begin
        Y:=2 / X;
        BBB:=Y * (Q5 + Y * (Q6 + Y * Q7));
        AAA:=Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + BBB)));
        BESSK0:=(Exp(-X) / Sqrt(X)) * AAA;
      end;
end;

Function BESSK(N:integer; X:real):real;
var
    TOX,BKM,BK,BKP:real; J:integer;
begin
    If N < 2 Then
        ShowMessage('bad argument N in BESSK');
    TOX:=2 / X;
    BKM:=BESSK0(X);
    BK:=BESSK1(X);
    For J:=1 To N - 1 do
    begin
        BKP:=BKM + J * TOX * BK;
        BKM:=BK;
        BK:=BKP;
    end;
    BESSK:=BK;
end;

end.