d4r28.frm

VB常用数值算法集2

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   5775
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   6150
   LinkTopic       =   "Form1"
   ScaleHeight     =   5775
   ScaleWidth      =   6150
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   4080
      TabIndex        =   0
      Top             =   4920
      Width           =   1455
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Private Sub Command1_ClicK()
    'PROGRAM D4R28
    'Driver for routine BESIAN
    Dim N(18), VALUE(18), F(18), A(18), X(18)
    Open "D:\VB常用数值算法集\DATA\FNCVAL.DAT" For Input As #1
    Do
      Line Input #1, Text$
    Loop While Text$ <> "Modified Bessel Function Ia+n"
    Line Input #1, NVAL$
    NVAL = Val(NVAL$)
    Print Text$
    Print
    Print "   N     A     IH        X          Actual       BESIAN"
    For I = 0 To NVAL - 1
        Line Input #1, DUM$
        N(I) = Val(Mid(DUM$, 1, 7))
        A(I) = Val(Mid(DUM$, 8, 8))
        X(I) = Val(Mid(DUM$, 16, 8))
        VALUE(I) = Val(Mid(DUM$, 24))
    Next I
    IH = 1
    Call BESIAN(X(0), A(0), 3, IH, F())
    For I = 0 To 3
        Print Tab(4); Format$(N(I), "#0");
        Print Tab(9); Format$(A(I), "0.#");
        Print Tab(16); Format$(IH, "#");
        Print Tab(22); Format$(X(I), "0.000000");
        Print Tab(34); Format$(VALUE(I), ".000000E+00");
        Print Tab(48); Format$(F(I), ".000000E+00");
    Next I
    Call BESIAN(X(4), A(4), 3, IH, F())
    For I = 4 To 7
        Print Tab(4); Format$(N(I), "#0");
        Print Tab(9); Format$(A(I), "0.#");
        Print Tab(16); Format$(IH, "#");
        Print Tab(22); Format$(X(I), "0.000000");
        Print Tab(34); Format$(VALUE(I), ".000000E+00");
        Print Tab(48); Format$(F(I - 4), ".000000E+00");
    Next I
    IH = -1
    Call BESIAN(X(8), A(8), 3, IH, F())
    For I = 8 To 10
        Print Tab(4); Format$(N(I), "#0");
        Print Tab(9); Format$(A(I), "0.#");
        Print Tab(16); Format$(IH, "#");
        Print Tab(22); Format$(X(I), "0.000000");
        Print Tab(34); Format$(VALUE(I), ".000000E+00");
        Print Tab(48); Format$(F(I - 7), ".000000E+00");
    Next I
    Call BESIAN(X(11), A(11), 3, IH, F())
    For I = 11 To 13
        Print Tab(4); Format$(N(I), "#0");
        Print Tab(9); Format$(A(I), "0.#");
        Print Tab(16); Format$(IH, "#");
        Print Tab(22); Format$(X(I), "0.000000");
        Print Tab(34); Format$(VALUE(I), ".000000E+00");
        Print Tab(48); Format$(F(I - 10), ".000000E+00");
    Next I
    Close #1
End Sub
Sub BESIAN(X, A, NM, IH, F())
    If IH = 1 Then
        Call IAP(X, A, NM, F())
        Exit Sub
    Else
        Call IAP(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 IAP(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#)) * Exp(X) / Sum
    D1 = 2.3026 * E + 1.3863
    If NMAX > 0 Then
        R = T(0.5 * D1 / NMAX) * NMAX
    Else
        R = 0#
    End If
    If X < D1 Then
        S = T(0.73576 * (D1 - X) / X) * 1.3591 * X
    Else
        S = 1.3591 * X
    End If
    If R <= S Then
        NU = 1 + Int(S)
    Else
        NU = 1 + Int(R)
    End If
1   N = 0
    AL = 1#
2   N = N + 1
    AL = AL * (N + 2# * A) / (N + 1#)
    If N < NU Then GoTo 2
    R = 0#
    S = 0#
3   R = 1# / (2# * (A + N) / X + R)
    AL = AL * (N + 1#) / (N + 2# * A)
    AM = 2# * (N + A) * AL
    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
   Function GAMMLN(XX)
      Dim COF(6)
      COF(1) = 76.18009173
      COF(2) = -86.50532033
      COF(3) = 24.01409822
      COF(4) = -1.231739516
      COF(5) = 0.00120858003
      COF(6) = -0.00000536382
      STP = 2.50662827465
      HALF = 0.5
      ONE = 1#
      FPF = 5.5
      X = XX - ONE
      TMP = X + FPF
      TMP = (X + HALF) * Log(TMP) - TMP
      SER = ONE
      For J = 1 To 6
          X = X + ONE
          SER = SER + COF(J) / X
    Next J
      GAMMLN = TMP + Log(STP * SER)
   End Function