Function T(Y:real):real;
var
AAA,Z,P:REAL;
begin
If Y 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 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 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;
