unit Unit2;

interface
Function GAMMLN(xx:real):real;
Function T(Y:real):real;
procedure JAP(X, A:real; NMAX:integer; VAR F:array of real);
procedure BESJAN(X, A:real; NM, IH:integer; VAR FF:array of real);


implementation
Function GAMMLN(xx:real):real;
CONST
  STP=2.50662827465; HALF=0.5; ONE=1.0; FPF=5.5;
var
  x,tmp,ser:double;
  j:integer;
  cof:array[1..6] of double;
BEGIN
  COF[1]:=76.18009173;      COF[2]:=-86.50532033;
  COF[3]:=24.01409822;      COF[4]:=-1.231739516;
  COF[5]:= 0.120858003e-2;  COF[6]:=-0.536382e-5;
  X:=XX-ONE;
  TMP:=X+FPF;
  TMP:=(X+HALF)*Ln(TMP)-TMP;
  SER:=ONE;
  For J:=1 To 6 do
    BEGIN
      X:=X+ONE;
      SER:=SER+COF[J]/X
    END;
  GAMMLN:=TMP+Ln(STP*SER);
end;

Function T(Y:real):real;
var
    AAA,Z,P:REAL;
begin
    If Y <= 10 Then
    begin
        AAA:=0.000057941 * Y - 0.00176148 * Y + 0.0208645;
        T:=((AAA * Y - 0.129013) * Y + 0.85777) * Y + 1.0125;
    end
    Else
    begin
        Z:=Ln(Y) - 0.775;
        P:=(0.775 - Ln(Z)) / (1 + Z);
        T:=Y / (P + 1) / Z;
    end;
end;

procedure JAP(X, A:real; NMAX:integer; VAR F:array of real);
label 1,2,3;
const
    MMAX = 10;   E = 4;
var
    N,NU,M,LI:Integer; EPS,SUM,D1,R,S,AM,AL:Real;
    FA,RR:array[0..10] of real;
begin
    EPS:= 0.5 * EXP(-E * LN(10));
    For N:=0 To NMAX do  FA[N]:=0;
    Sum:=Exp(GAMMLN(1 + A));
    Sum:=Exp(A * Ln(X / 2)) / Sum;
    D1:=2.3026 * E + 1.3863;
    If NMAX > 0 Then
        R:=T(0.5 * D1 / NMAX) * NMAX
    Else
        R:=0;
    S:=T(0.73576 * D1 / X) * 1.3591 * X;
    If R <= S Then
        NU:=1 + Trunc(S)
    Else
        NU:=1 + Trunc(R);
1:  M:=0;
    AL:=1;
    LI:=Trunc(NU / 2);
2:  M:=M + 1;
    AL:=AL * (M + A) / (M + 1);
    If M < LI Then GoTo 2;
    N:=2 * M;
    R:=0;
    S:=0;
3:  R:=1 / (2 * (A + N) / X - R);
    If Trunc(N / 2) <> N / 2 Then
        AM:=0
    Else
    begin
        AL:=AL * (N + 2) / (N + 2 * A);
        AM:=AL * (N + A);
    end;
    S:=R * (AM + S);
    If N <= NMAX Then RR[N - 1]:=R;
    N:=N - 1;
    If N >= 1 Then GoTo 3;
    F[0]:=Sum / (1 + S);
    For N:=0 To NMAX - 1 do F[N + 1]:=RR[N] * F[N];
    For N:=0 To NMAX do
    begin
        If Abs((F[N] - FA[N]) / F[N]) > EPS Then
        begin
            For M:=0 To NMAX do FA[M]:=F[M];
            NU:=NU + 5;
            GoTo 1;
        end;
    end;
end;

procedure BESJAN(X, A:real; NM, IH:integer; var FF:array of real);
var
    I:integer;
begin
    If IH =1 Then
      JAP(X, A, NM, FF)
    Else
    begin
        JAP(X, A, 1, FF);
        FF[1]:= 2 * A * FF[0] / X - FF[1];
        For I:=1 To NM - 1 do
          FF[I + 1]:= 2 * (A - I) * FF[I] / X - FF[I - 1];
    end;
end;

end.