Sub BESJAN(X, A, NM, IH, F())
If IH = 1 Then
Call JAP(X, A, NM, F())
Exit Sub
Else
Call JAP(X, A, 1, F())
F(1) = 2# * A * F(0) / X - F(1)
For I = 1 To NM - 1
F(I + 1) = 2# * (A - I) * F(I) / X - F(I - 1)
Next I
Exit Sub
End If
End Sub
Sub JAP(X, A, NMAX, F())
MMAX = 10
E = 4
EPS = 0.5 * 10 ^ (-E)
Dim FA(10), RR(10)
For N = 0 To NMAX
FA(N) = 0#
Next N
Sum = Exp(GAMMLN(1# + A))
Sum = Exp(A * Log(X / 2#)) / Sum
D1 = 2.3026 * E + 1.3863
If NMAX > 0 Then
R = T(0.5 * D1 / NMAX) * NMAX
Else
R = 0#
End If
S = T(0.73576 * D1 / X) * 1.3591 * X
If R NU = 1 + Int(S)
Else
NU = 1 + Int(R)
End If
1 M = 0
AL = 1#
LI = Int(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 Int(N / 2) N / 2# Then
AM = 0#
Else
AL = AL * (N + 2) / (N + 2# * A)
AM = AL * (N + A)
End If
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
F(N + 1) = RR(N) * F(N)
Next N
For N = 0 To NMAX
If Abs((F(N) - FA(N)) / F(N)) > EPS Then
For M = 0 To NMAX
FA(M) = F(M)
Next M
NU = NU + 5
GoTo 1
End If
Next N
End Sub
Function T(Y)
If Y AAA = 0.000057941 * Y - 0.00176148 * Y + 0.0208645
T = ((AAA * Y - 0.129013) * Y + 0.85777) * Y + 1.0125
Else
Z = Log(Y) - 0.775
P = (0.775 - Log(Z)) / (1# + Z)
T = Y / (P + 1#) / Z
End If
End Function
