数学算法

				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;