d9r8b.frm

数据拟合源程序

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   6315
   ClientLeft      =   1455
   ClientTop       =   945
   ClientWidth     =   7410
   LinkTopic       =   "Form1"
   ScaleHeight     =   6315
   ScaleWidth      =   7410
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   5880
      TabIndex        =   0
      Top             =   5640
      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 D9R8B
    'Driver for routine MRQCOF
    NPT = 100
    MA = 6
    SPREAD = 0.1
    Dim X(100), Y(100), SIG(100), A(6), LISTA(6)
    Dim COVAR(6, 6), ALPHA(6, 6), GUES(6), BETA(6)
    A(1) = 5#: A(2) = 2#: A(3) = 3#: A(4) = 2#: A(5) = 5#: A(6) = 3#
    GUES(1) = 4.9: GUES(2) = 2.1: GUES(3) = 2.9
    GUES(4) = 2.1: GUES(5) = 4.9: GUES(6) = 3.1
    IDUM& = -911
    'First try a sum of two Gaussians
    For I = 1 To 100
        X(I) = 0.1 * I
        Y(I) = 0#
        For J = 1 To 4 Step 3
            Y(I) = Y(I) + A(J) * Exp(-((X(I) - A(J + 1)) / A(J + 2)) ^ 2)
        Next J
        Y(I) = Y(I) * (1# + SPREAD * GASDEV(IDUM&))
        SIG(I) = SPREAD * Y(I)
    Next I
    MFIT = MA
    For I = 1 To MFIT
        LISTA(I) = I
    Next I
    For I = 1 To MA
        A(I) = GUES(I)
    Next I
    Call MRQCOF(X(), Y(), SIG(), NPT, A(), MA, LISTA(), MFIT, ALPHA(), BETA(), MA, CHISQ)
    Print Tab(5)
    Print Tab(5); "matrix alpha"
    Print Tab(5)
    For I = 1 To MA
        For J = 1 To MA
            Print Tab((J - 1) * 12); Format$(ALPHA(I, J), "####.###0");
        Next J
        Print
    Next I
    Print Tab(5)
    Print Tab(5); "vector beta"
    Print Tab(5)
    For I = 1 To MA
        Print Format$(BETA(I), "###.###0"),
    Next I
    Print Tab(5)
    Print Tab(5); "Chi-squared: "; Format$(CHISQ, "###.####")
    Print Tab(5)
    'Next fix one line and improve the other
    For I = 1 To 3
        LISTA(I) = I + 3
    Next I
    MFIT = 3
    For I = 1 To MA
        A(I) = GUES(I)
    Next I
    Call MRQCOF(X(), Y(), SIG(), NPT, A(), MA, LISTA(), MFIT, ALPHA(), BETA(), MA, CHISQ)
    Print Tab(5); "matrix alpha"
    Print Tab(5)
    For I = 1 To MFIT
        For J = 1 To MFIT
            Print Tab((J - 1) * 12); Format$(ALPHA(I, J), "####.###0");
        Next J
        Print
    Next I
    Print Tab(5)
    Print Tab(5); "vector beta"
    Print Tab(5)
    For I = 1 To MFIT
        Print Format$(BETA(I), "###.###0"),
    Next I
    Print Tab(5)
    Print Tab(5); "Chi-squared: "; Format$(CHISQ, "###.###0")
End Sub
Sub MRQCOF(X(), Y(), SIG(), NDATA, A(), MA, LISTA(), MFIT, ALPHA(), BETA(), NALP, CHISQ)
    Dim DYDA(20)
    For J = 1 To MFIT
        For K = 1 To J
            ALPHA(J, K) = 0#
        Next K
        BETA(J) = 0#
    Next J
    CHISQ = 0#
    For I = 1 To NDATA
        Call FGAUSS(X(I), A(), YMOD, DYDA(), MA)
        SIG2I = 1# / (SIG(I) * SIG(I))
        DY = Y(I) - YMOD
        For J = 1 To MFIT
            WT = DYDA(LISTA(J)) * SIG2I
            For K = 1 To J
                ALPHA(J, K) = ALPHA(J, K) + WT * DYDA(LISTA(K))
            Next K
            BETA(J) = BETA(J) + DY * WT
        Next J
        CHISQ = CHISQ + DY * DY * SIG2I
    Next I
    For J = 2 To MFIT
        For K = 1 To J - 1
            ALPHA(K, J) = ALPHA(J, K)
        Next K
    Next J
End Sub
Sub FGAUSS(X, A(), Y, DYDA(), NA)
    Y = 0#
    For I = 1 To NA - 1 Step 3
        ARG = (X - A(I + 1)) / A(I + 2)
        EX = Exp(-ARG ^ 2)
        FAC = A(I) * EX * 2# * ARG
        Y = Y + A(I) * EX
        DYDA(I) = EX
        DYDA(I + 1) = FAC / A(I + 2)
        DYDA(I + 2) = FAC * ARG / A(I + 2)
    Next I
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