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 <= S Then
        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 <= 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
        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 <= 10# Then
        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