?? d9r2.frm
字號(hào):
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
?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -