📄 d9r4.frm
字号:
VERSION 5.00
Begin VB.Form Form1
Caption = "Form1"
ClientHeight = 4920
ClientLeft = 60
ClientTop = 345
ClientWidth = 4680
LinkTopic = "Form1"
ScaleHeight = 4920
ScaleWidth = 4680
StartUpPosition = 3 'Windows Default
Begin VB.CommandButton Command1
Caption = "Command1"
Height = 375
Left = 3000
TabIndex = 0
Top = 4320
Width = 1335
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 D9R4
'Driver for routine SVDFIT
NPT = 100
SPREAD = 0.02
NPOL = 5
Dim X(100), Y(100), SIG(100), A(5), CVM(5, 5)
Dim U(100, 5), V(5, 5), W(5)
'Polynomial fit
IDUM& = -911
MP = NPT
NP = NPOL
For I = 1 To NPT
X(I) = 0.02 * I
Y(I) = 1# + X(I) * (2# + X(I) * (3# + X(I) * (4# + X(I) * 5#)))
Y(I) = Y(I) + SPREAD * GASDEV(IDUM&)
SIG(I) = SPREAD
Next I
Call SVDFIT(X(), Y(), SIG(), NPT, A(), NPOL, U(), V(), W(), MP, NP, CHISQ, "FPOLY")
Call SVDVAR(V(), NPOL, NP, W(), CVM(), NPOL)
Print
Print Tab(5); "Polynomial fit"
Print
For I = 1 To NPOL
Print Tab(5); Format$(A(I), "#.##0000"), " +-",
Print Format$(Sqr(CVM(I, I)), ".#####0")
Next I
Print
Print Tab(5); "Chi-squared "; Format$(CHISQ, "###.####00")
Call SVDFIT(X(), Y(), SIG(), NPT, A(), NPOL, U(), V(), W(), MP, NP, CHISQ, "FLEG")
Call SVDVAR(V(), NPOL, NP, W(), CVM(), NPOL)
Print
Print Tab(5); "Legendre polynomial fit"
For I = 1 To NPOL
Print Tab(5); Format$(A(I), "#.##0000"), " +-",
Print Format$(Sqr(CVM(I, I)), ".#####0")
Next I
Print
Print Tab(5); "Chi-squared "; Format$(CHISQ, "###.####00")
End Sub
Sub SVDFIT(X(), Y(), SIG(), NDATA, A(), MA, U(), V(), W(), MP, NP, CHISQ, FUNCS$)
TOL = 0.00001
Dim B(1000), AFUNC(50)
For I = 1 To NDATA
If FUNCS$ = "FPOLY" Then Call FPOLY(X(I), AFUNC(), MA)
If FUNCS$ = "FLEG" Then Call FLEG(X(I), AFUNC(), MA)
TMP = 1# / SIG(I)
For J = 1 To MA
U(I, J) = AFUNC(J) * TMP
Next J
B(I) = Y(I) * TMP
Next I
Call SVDCMP(U(), NDATA, MA, W(), V())
WMAX = 0#
For J = 1 To MA
If W(J) > WMAX Then WMAX = W(J)
Next J
THRESH = TOL * WMAX
For J = 1 To MA
If W(J) < THRESH Then W(J) = 0#
Next J
Call SVBKSB(U(), W(), V(), NDATA, MA, B(), A())
CHISQ = 0#
For I = 1 To NDATA
If FUNCS$ = "FPOLY" Then Call FPOLY(X(I), AFUNC(), MA)
If FUNCS$ = "FLEG" Then Call FLEG(X(I), AFUNC(), MA)
SUM1 = 0#
For J = 1 To MA
SUM1 = SUM1 + A(J) * AFUNC(J)
Next J
CHISQ = CHISQ + ((Y(I) - SUM1) / SIG(I)) ^ 2
Next I
End Sub
Sub FPOLY(X, P(), NP)
P(1) = 1#
For J = 2 To NP
P(J) = P(J - 1) * X
Next J
End Sub
Sub FLEG(X, PL(), NL)
PL(1) = 1#
PL(2) = X
If NL > 2 Then
TWOX = 2# * X
F2 = X
D = 1#
For J = 3 To NL
F1 = D
F2 = F2 + TWOX
D = D + 1#
PL(J) = (F2 * PL(J - 1) - F1 * PL(J - 2)) / D
Next J
End If
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 SVDVAR(V(), MA, NP, W(), CVM(), NCVM)
Dim WTI(20)
For I = 1 To MA
WTI(I) = 0#
If W(I) <> 0# Then WTI(I) = 1# / (W(I) * W(I))
Next I
For I = 1 To MA
For J = 1 To I
SUM1 = 0#
For K = 1 To MA
SUM1 = SUM1 + V(I, K) * V(J, K) * WTI(K)
Next K
CVM(I, J) = SUM1
CVM(J, I) = SUM1
Next J
Next I
End Sub
Sub SVDCMP(A(), M, N, W(), V())
Dim RV1(100)
If M < N Then Print "You must augment A with extra zero rows."
G = 0#
SCALE1 = 0#
ANORM = 0#
For I = 1 To N
L = I + 1
RV1(I) = SCALE1 * G
G = 0#
S = 0#
SCALE1 = 0#
If I <= M Then
For K = I To M
SCALE1 = SCALE1 + Abs(A(K, I))
Next K
If SCALE1 <> 0# Then
For K = I To M
A(K, I) = A(K, I) / SCALE1
S = S + A(K, I) * A(K, I)
Next K
F = A(I, I)
G = -Sqr(S) * Sgn(F)
H = F * G - S
A(I, I) = F - G
If I <> N Then
For J = L To N
S = 0#
For K = I To M
S = S + A(K, I) * A(K, J)
Next K
F = S / H
For K = I To M
A(K, J) = A(K, J) + F * A(K, I)
Next K
Next J
End If
For K = I To M
A(K, I) = SCALE1 * A(K, I)
Next K
End If
End If
W(I) = SCALE1 * G
G = 0#
S = 0#
SCALE1 = 0#
If I <= M And I <> N Then
For K = L To N
SCALE1 = SCALE1 + Abs(A(I, K))
Next K
If SCALE1 <> 0# Then
For K = L To N
A(I, K) = A(I, K) / SCALE1
S = S + A(I, K) * A(I, K)
Next K
F = A(I, L)
G = -Sqr(S) * Sgn(F)
H = F * G - S
A(I, L) = F - G
For K = L To N
RV1(K) = A(I, K) / H
Next K
If I <> M Then
For J = L To M
S = 0#
For K = L To N
S = S + A(J, K) * A(I, K)
Next K
For K = L To N
A(J, K) = A(J, K) + S * RV1(K)
Next K
Next J
End If
For K = L To N
A(I, K) = SCALE1 * A(I, K)
Next K
End If
End If
If ANORM > Abs(W(I)) + Abs(RV1(I)) Then
ANORM = ANORM
Else
ANORM = Abs(W(I)) + Abs(RV1(I))
End If
Next I
For I = N To 1 Step -1
If I < N Then
If G <> 0# Then
For J = L To N
V(J, I) = (A(I, J) / A(I, L)) / G
Next J
For J = L To N
S = 0#
For K = L To N
S = S + A(I, K) * V(K, J)
Next K
For K = L To N
V(K, J) = V(K, J) + S * V(K, I)
Next K
Next J
End If
For J = L To N
V(I, J) = 0#
V(J, I) = 0#
Next J
End If
V(I, I) = 1#
G = RV1(I)
L = I
Next I
For I = N To 1 Step -1
L = I + 1
G = W(I)
If I < N Then
For J = L To N
A(I, J) = 0#
Next J
End If
If G <> 0# Then
G = 1# / G
If I <> N Then
For J = L To N
S = 0#
For K = L To M
S = S + A(K, I) * A(K, J)
Next K
F = (S / A(I, I)) * G
For K = I To M
A(K, J) = A(K, J) + F * A(K, I)
Next K
Next J
End If
For J = I To M
A(J, I) = A(J, I) * G
Next J
Else
For J = I To M
A(J, I) = 0#
Next J
End If
A(I, I) = A(I, I) + 1#
Next I
For K = N To 1 Step -1
For ITS = 1 To 30
For L = K To 1 Step -1
NM = L - 1
If Abs(RV1(L)) + ANORM = ANORM Then GoTo 2
If Abs(W(NM)) + ANORM = ANORM Then GoTo 1
Next L
1 C = 0#
S = 1#
For I = L To K
F = S * RV1(I)
If Abs(F) + ANORM <> ANORM Then
G = W(I)
H = Sqr(F * F + G * G)
W(I) = H
H = 1# / H
C = (G * H)
S = -(F * H)
For J = 1 To M
Y = A(J, NM)
Z = A(J, I)
A(J, NM) = (Y * C) + (Z * S)
A(J, I) = -(Y * S) + (Z * C)
Next J
End If
Next I
2 Z = W(K)
If L = K Then
If Z < 0# Then
W(K) = -Z
For J = 1 To N
V(J, K) = -V(J, K)
Next J
End If
GoTo 3
End If
If ITS = 30 Then Print "No convergence in 30 iterations"
X = W(L)
NM = K - 1
Y = W(NM)
G = RV1(NM)
H = RV1(K)
F = ((Y - Z) * (Y + Z) + (G - H) * (G + H)) / (2# * H * Y)
G = Sqr(F * F + 1#)
F = ((X - Z) * (X + Z) + H * ((Y / (F + Abs(G) * Sgn(F))) - H)) / X
C = 1#
S = 1#
For J = L To NM
I = J + 1
G = RV1(I)
Y = W(I)
H = S * G
G = G * C
Z = Sqr(F * F + H * H)
RV1(J) = Z
C = F / Z
S = H / Z
F = (X * C) + (G * S)
G = -(X * S) + (G * C)
H = Y * S
Y = Y * C
For NM = 1 To N
X = V(NM, J)
Z = V(NM, I)
V(NM, J) = (X * C) + (Z * S)
V(NM, I) = -(X * S) + (Z * C)
Next NM
Z = Sqr(F * F + H * H)
W(J) = Z
If Z <> 0# Then
Z = 1# / Z
C = F * Z
S = H * Z
End If
F = (C * G) + (S * Y)
X = -(S * G) + (C * Y)
For NM = 1 To M
Y = A(NM, J)
Z = A(NM, I)
A(NM, J) = (Y * C) + (Z * S)
A(NM, I) = -(Y * S) + (Z * C)
Next NM
Next J
RV1(L) = 0#
RV1(K) = F
W(K) = X
Next ITS
3 AAAAA = 1
Next K
End Sub
Sub SVBKSB(U(), W(), V(), M, N, B(), X())
Dim TMP(100)
For J = 1 To N
S = 0#
If W(J) <> 0# Then
For I = 1 To M
S = S + U(I, J) * B(I)
Next I
S = S / W(J)
End If
TMP(J) = S
Next J
For J = 1 To N
S = 0#
For JJ = 1 To N
S = S + V(J, JJ) * TMP(JJ)
Next JJ
X(J) = S
Next J
End Sub
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -