d9r2.frm

常用的数值算法的VB程序

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   7050
   ClientLeft      =   2265
   ClientTop       =   345
   ClientWidth     =   5130
   LinkTopic       =   "Form1"
   ScaleHeight     =   7050
   ScaleWidth      =   5130
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3720
      TabIndex        =   0
      Top             =   6480
      Width           =   1215
   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 D9R2
    'Driver for routine LFIT
    NPT = 100
    SPREAD = 0.1
    NTERM = 3
    Dim X(100), Y(100), SIG(100), A(3), LISTA(3), COVAR(3, 3)
    IDUM& = -911
    For I = 1 To NPT
        X(I) = 0.1 * I
        Y(I) = NTERM
        For J = NTERM - 1 To 1 Step -1
            Y(I) = J + Y(I) * X(I)
        Next J
        Y(I) = Y(I) + SPREAD * GASDEV(IDUM&)
        SIG(I) = SPREAD
    Next I
    MFIT = NTERM
    For I = 1 To MFIT
        LISTA(I) = I
    Next I
    Call LFIT(X(), Y(), SIG(), NPT, A(), NTERM, LISTA(), MFIT, COVAR(), NTERM, CHISQ)
    Print Tab(8); "PARAMETER           Uncertainty"
    For I = 1 To NTERM
        Print Tab(5); "A("; I; ")= "; Format$(A(I), "##.000000"),
        Print Format$(Sqr(COVAR(I, I)), "##.000000")
    Next I
    Print Tab(5)
    Print Tab(5); "Chi-squared =  "; Format$(CHISQ, ".###000E+00")
    Print Tab(5)
    Print Tab(5); "Full covariance matrix"
    For I = 1 To NTERM
        For J = 1 To NTERM
            Print Tab(8 + (J - 1) * 13); Format$(COVAR(I, J), ".00E+00");
        Next J
    Next I
    'Now test the LISTA feature
    For I = 1 To NTERM
        LISTA(I) = NTERM + 1 - I
    Next I
    Call LFIT(X(), Y(), SIG(), NPT, A(), NTERM, LISTA(), MFIT, COVAR(), NTERM, CHISQ)
    Print Tab(5)
    Print Tab(8); "PARAMETER           Uncertainty"
    For I = 1 To NTERM
        Print Tab(5); "A("; I; ")= "; Format$(A(I), "##.000000"),
        Print Format$(Sqr(COVAR(I, I)), "##.000000")
    Next I
    Print Tab(5)
    Print Tab(5); "Chi-squared =  "; Format$(CHISQ, ".###000E+00")
    Print Tab(5)
    Print Tab(5); "Full covariance matrix"
    For I = 1 To NTERM
        For J = 1 To NTERM
            Print Tab(8 + (J - 1) * 13); Format$(COVAR(I, J), ".00E+00");
        Next J
    Next I
    'Now check results of restricting fit parameters
    II = 1
    For I = 1 To NTERM
        AAA = I - Int(I / 2) * 2
        If AAA = 1 Then
            LISTA(II) = I
            II = II + 1
        End If
    Next
    MFIT = II - 1
    Call LFIT(X(), Y(), SIG(), NPT, A(), NTERM, LISTA(), MFIT, COVAR(), NTERM, CHISQ)
    Print Tab(5)
    Print Tab(8); "PARAMETER           Uncertainty"
    For I = 1 To NTERM
        Print Tab(5); "A("; I; ")= "; Format$(A(I), "##.000000"),
        Print Format$(Sqr(COVAR(I, I)), "##.000000")
    Next
    Print Tab(5)
    Print Tab(5); "Chi-squared = "; Format$(CHISQ, ".###000E+00")
    Print Tab(5)
    Print Tab(5); "Full covariance matrix"
    For I = 1 To NTERM
        For J = 1 To NTERM
            Print Tab(8 + (J - 1) * 13); Format$(COVAR(I, J), ".00E+00");
        Next J
    Next I
End Sub
Sub FUNCS(X, AFUNC, MA)
    AFUNC(1) = 1#
    For I = 2 To MA
        AFUNC(I) = X * AFUNC(I - 1)
    Next I
End Sub
Sub LFIT(X(), Y(), SIG(), NDATA, A(), MA, LISTA(), MFIT, COVAR(), NCVM, CHISQ)
    Dim BETA(50), AFUNC(50)
    KK = MFIT + 1
    For J = 1 To MA
        IHIT = 0
        For K = 1 To MFIT
            If LISTA(K) = J Then IHIT = IHIT + 1
        Next K
        If IHIT = 0 Then
            LISTA(KK) = J
            KK = KK + 1
        ElseIf IHIT > 1 Then
            Print " Improper set in LISTA"
        End If
    Next J
    If KK <> (MA + 1) Then Print " Improper set in LISTA"
    For J = 1 To MFIT
        For K = 1 To MFIT
            COVAR(J, K) = 0#
        Next K
        BETA(J) = 0#
    Next J
    For I = 1 To NDATA
        Call FUNCS(X(I), AFUNC, MA)
        YM = Y(I)
        If MFIT < MA Then
            For J = MFIT + 1 To MA
                YM = YM - A(LISTA(J)) * AFUNC(LISTA(J))
            Next J
        End If
        SIG2I = 1# / SIG(I) ^ 2
        For J = 1 To MFIT
            WT = AFUNC(LISTA(J)) * SIG2I
            For K = 1 To J
                COVAR(J, K) = COVAR(J, K) + WT * AFUNC(LISTA(K))
           Next K
            BETA(J) = BETA(J) + YM * WT
        Next J
    Next I
    If MFIT > 1 Then
        For J = 2 To MFIT
            For K = 1 To J - 1
                COVAR(K, J) = COVAR(J, K)
            Next K
        Next J
    End If
    Call GAUSSJ(COVAR(), MFIT, BETA())
    For J = 1 To MFIT
        A(LISTA(J)) = BETA(J)
    Next J
    CHISQ = 0#
    For I = 1 To NDATA
        Call FUNCS(X(I), AFUNC, MA)
        Sum = 0#
        For J = 1 To MA
            Sum = Sum + A(J) * AFUNC(J)
        Next J
        CHISQ = CHISQ + ((Y(I) - Sum) / SIG(I)) ^ 2
    Next I
    Call COVSRT(COVAR, NCVM, MA, LISTA, MFIT)
End Sub
Function GASDEV(IDUM&)
    Static ISET, GSET
    If ISET = 0 Then
      Do
        V1 = 2# * RAN1(IDUM&) - 1#
        V2 = 2# * RAN1(IDUM&) - 1#
        R = V1 ^ 2 + V2 ^ 2
      Loop While R >= 1# Or R = 0
      FAC = Sqr(-2# * Log(R) / R)
      GSET = V1 * FAC
      GASDEV = V2 * FAC
      ISET = 1
    Else
      GASDEV = GSET
      ISET = 0
    End If
End Function
Static Function RAN1(IDUM&)
    Dim R(97)
    M1& = 259200: IA1& = 7141: IC1& = 54773: RM1 = 0.0000038580247
    M2& = 134456: IA2& = 8121: IC2& = 28411: RM2 = 0.0000074373773
    M3& = 243000: IA3& = 4561: IC3& = 51349
    If IDUM& < 0 Or IFF = 0 Then
        IFF = 1
        IX1& = (IC1& - IDUM&) Mod M1&
        IX1& = (IA1& * IX1& + IC1&) Mod M1&
        IX2& = IX1& Mod M2&
        IX1& = (IA1& * IX1& + IC1&) Mod M1&
        IX3& = IX1& Mod M3&
        For J = 1 To 97
            IX1& = (IA1& * IX1& + IC1&) Mod M1&
            IX2& = (IA2& * IX2& + IC2&) Mod M2&
            R(J) = (CSng(IX1&) + CSng(IX2&) * RM2) * RM1
        Next J
        IDUM& = 1
    End If
    IX1& = (IA1& * IX1& + IC1&) Mod M1&
    IX2& = (IA2& * IX2& + IC2&) Mod M2&
    IX3& = (IA3& * IX3& + IC3&) Mod M3&
    J = 1 + Int((97 * IX3&) / M3&)
    If J > 97 Or J < 1 Then Print "Abnormal exit": Exit Function
    RAN1 = R(J)
    R(J) = (CSng(IX1&) + CSng(IX2&) * RM2) * RM1
End Function
Sub COVSRT(COVAR(), NCVM, MA, LISTA(), MFIT)
      For J = 1 To MA - 1
          For I = J + 1 To MA
              COVAR(I, J) = 0#
          Next I
      Next J
      For I = 1 To MFIT - 1
          For J = I + 1 To MFIT
              If LISTA(J) > LISTA(I) Then
                  COVAR(LISTA(J), LISTA(I)) = COVAR(I, J)
              Else
                  COVAR(LISTA(I), LISTA(J)) = COVAR(I, J)
              End If
          Next J
      Next I
      SWAP = COVAR(1, 1)
      For J = 1 To MA
          COVAR(1, J) = COVAR(J, J)
          COVAR(J, J) = 0#
      Next J
      COVAR(LISTA(1), LISTA(1)) = SWAP
      For J = 2 To MFIT
          COVAR(LISTA(J), LISTA(J)) = COVAR(1, J)
      Next J
      For J = 2 To MA
          For I = 1 To J - 1
              COVAR(I, J) = COVAR(J, I)
          Next I
      Next J
   End Sub
Sub GAUSSJ(A(), N, B())
      Dim IPIV(50), INDXR(50), INDXC(50)
      For J = 1 To N
          IPIV(J) = 0
      Next J
      For I = 1 To N
          BIG = 0#
          For J = 1 To N
              If IPIV(J) <> 1 Then
                  For K = 1 To N
                    If IPIV(K) = 0 Then
                      If Abs(A(J, K)) >= BIG Then
                          BIG = Abs(A(J, K))
                          IROW = J
                          ICOL = K
                      End If
                    ElseIf IPIV(K) > 1 Then
                      Print "Singular matrix"
                    End If
                  Next K
              End If
          Next J
          IPIV(ICOL) = IPIV(ICOL) + 1
          If IROW <> ICOL Then
              For L = 1 To N
                  DUM = A(IROW, L)
                  A(IROW, L) = A(ICOL, L)
                  A(ICOL, L) = DUM
              Next L
              DUM = B(IROW)
              B(IROW) = B(ICOL)
              B(ICOL) = DUM
          End If
          INDXR(I) = IROW
          INDXC(I) = ICOL
          If A(ICOL, ICOL) = 0# Then Print "Singular matrix."
          PIVINV = 1# / A(ICOL, ICOL)
          A(ICOL, ICOL) = 1#
          For L = 1 To N
              A(ICOL, L) = A(ICOL, L) * PIVINV
          Next L
          B(ICOL) = B(ICOL) * PIVINV
          For LL = 1 To N
              If LL <> ICOL Then
                  DUM = A(LL, ICOL)
                  A(LL, ICOL) = 0#
                  For L = 1 To N
                      A(LL, L) = A(LL, L) - A(ICOL, L) * DUM
                  Next L
                  B(LL) = B(LL) - B(ICOL) * DUM
              End If
          Next LL
      Next I
      For L = N To 1 Step -1
           If INDXR(L) <> INDXC(L) Then
               For K = 1 To N
                   DUM = A(K, INDXR(L))
                   A(K, INDXR(L)) = A(K, INDXC(L))
                   A(K, INDXC(L)) = DUM
               Next K
           End If
      Next L
End Sub