www.pudn.com > 603Grid.rar > FrmGrid.frm, change:2003-09-13,size:71965b
VERSION 5.00
Object = "{F9043C88-F6F2-101A-A3C9-08002B2F49FB}#1.2#0"; "COMDLG32.OCX"
Begin VB.Form FrmGrid
BorderStyle = 3 'Fixed Dialog
Caption = "网格化程序"
ClientHeight = 8268
ClientLeft = 48
ClientTop = 336
ClientWidth = 11412
BeginProperty Font
Name = "宋体"
Size = 11.4
Charset = 134
Weight = 400
Underline = 0 'False
Italic = 0 'False
Strikethrough = 0 'False
EndProperty
LinkTopic = "Form1"
MaxButton = 0 'False
MinButton = 0 'False
ScaleHeight = 689
ScaleMode = 3 'Pixel
ScaleWidth = 951
ShowInTaskbar = 0 'False
StartUpPosition = 2 '屏幕中心
Begin VB.PictureBox Picture1
Height = 8052
Left = 120
ScaleHeight = 667
ScaleMode = 3 'Pixel
ScaleWidth = 927
TabIndex = 0
Top = 120
Width = 11172
Begin VB.PictureBox PictureGrid
Appearance = 0 'Flat
AutoRedraw = -1 'True
BackColor = &H80000005&
Enabled = 0 'False
BeginProperty Font
Name = "MS Sans Serif"
Size = 7.8
Charset = 0
Weight = 400
Underline = 0 'False
Italic = 0 'False
Strikethrough = 0 'False
EndProperty
ForeColor = &H80000008&
Height = 7680
Left = 3240
ScaleHeight = 7656
ScaleWidth = 7656
TabIndex = 20
Top = 120
Width = 7680
End
Begin VB.PictureBox PictureTools
Appearance = 0 'Flat
BackColor = &H80000000&
BeginProperty Font
Name = "MS Sans Serif"
Size = 7.8
Charset = 0
Weight = 400
Underline = 0 'False
Italic = 0 'False
Strikethrough = 0 'False
EndProperty
ForeColor = &H80000008&
Height = 7692
Left = 120
ScaleHeight = 7668
ScaleWidth = 2988
TabIndex = 1
Top = 120
Width = 3012
Begin VB.ListBox ListGrid
Appearance = 0 'Flat
Enabled = 0 'False
Height = 1848
ItemData = "FrmGrid.frx":0000
Left = 120
List = "FrmGrid.frx":0019
TabIndex = 23
Top = 3720
Width = 2772
End
Begin VB.TextBox TextXStep
Appearance = 0 'Flat
Enabled = 0 'False
Height = 336
Left = 120
TabIndex = 12
Top = 2280
Width = 1332
End
Begin VB.TextBox TextYStep
Appearance = 0 'Flat
Enabled = 0 'False
Height = 336
Left = 1560
TabIndex = 11
Top = 2280
Width = 1332
End
Begin VB.TextBox TextXNX
Appearance = 0 'Flat
Enabled = 0 'False
Height = 336
Left = 120
TabIndex = 10
Top = 3000
Width = 1332
End
Begin VB.TextBox TextYNY
Appearance = 0 'Flat
Enabled = 0 'False
Height = 336
Left = 1560
TabIndex = 9
Top = 3000
Width = 1332
End
Begin VB.CommandButton CommandGridOK
Caption = "开始网格化"
Enabled = 0 'False
BeginProperty Font
Name = "宋体"
Size = 11.4
Charset = 134
Weight = 700
Underline = 0 'False
Italic = 0 'False
Strikethrough = 0 'False
EndProperty
Height = 375
Left = 120
TabIndex = 8
Top = 6000
Width = 2775
End
Begin VB.CommandButton CommandGridCacel
Caption = "退出"
BeginProperty Font
Name = "宋体"
Size = 11.4
Charset = 134
Weight = 700
Underline = 0 'False
Italic = 0 'False
Strikethrough = 0 'False
EndProperty
Height = 375
Left = 120
TabIndex = 7
Top = 6720
Width = 2775
End
Begin VB.TextBox TextXmin
Appearance = 0 'Flat
Enabled = 0 'False
Height = 336
Left = 120
TabIndex = 6
Top = 840
Width = 1335
End
Begin VB.TextBox TextXmax
Appearance = 0 'Flat
Enabled = 0 'False
Height = 336
Left = 120
TabIndex = 5
Top = 1560
Width = 1332
End
Begin VB.TextBox TextYmin
Appearance = 0 'Flat
Enabled = 0 'False
Height = 336
Left = 1560
TabIndex = 4
Top = 840
Width = 1335
End
Begin VB.TextBox TextYmax
Appearance = 0 'Flat
Enabled = 0 'False
Height = 336
Left = 1560
TabIndex = 3
Top = 1560
Width = 1332
End
Begin VB.CommandButton CommandContouMeshWang
Caption = "打开文件..."
BeginProperty Font
Name = "宋体"
Size = 11.4
Charset = 134
Weight = 700
Underline = 0 'False
Italic = 0 'False
Strikethrough = 0 'False
EndProperty
Height = 375
Left = 120
MousePointer = 1 'Arrow
TabIndex = 2
Top = 120
Width = 2772
End
Begin VB.Label LabelYmax
Alignment = 2 'Center
Caption = "Y最大值"
Enabled = 0 'False
Height = 252
Left = 1560
TabIndex = 22
Top = 1320
Width = 1332
End
Begin VB.Label LabelXmax
Alignment = 2 'Center
Caption = "X最大值"
Enabled = 0 'False
Height = 252
Left = 120
TabIndex = 21
Top = 1320
Width = 1332
End
Begin VB.Label LabelXStep
Alignment = 2 'Center
Caption = "X方向步长"
Enabled = 0 'False
Height = 252
Left = 120
TabIndex = 19
Top = 2040
Width = 1332
End
Begin VB.Label LabelYStep
Alignment = 2 'Center
Caption = "Y方向步长"
Enabled = 0 'False
Height = 252
Left = 1560
TabIndex = 18
Top = 2040
Width = 1332
End
Begin VB.Label LabelXNX
Alignment = 2 'Center
Caption = "X方向网格数"
Enabled = 0 'False
Height = 252
Left = 120
TabIndex = 17
Top = 2760
Width = 1332
End
Begin VB.Label LabelYNY
Alignment = 2 'Center
Caption = "Y方向网格数"
Enabled = 0 'False
Height = 252
Left = 1560
TabIndex = 16
Top = 2760
Width = 1332
End
Begin VB.Label LabelGrid
Alignment = 2 'Center
Caption = "网格化方法"
Height = 252
Left = 120
TabIndex = 15
Top = 3480
Width = 2772
End
Begin VB.Label LabelXmin
Alignment = 2 'Center
Caption = "X最小值"
Enabled = 0 'False
Height = 252
Left = 120
TabIndex = 14
Top = 600
Width = 1332
End
Begin VB.Label LabelYmin
Alignment = 2 'Center
Caption = "Y最小值"
Enabled = 0 'False
Height = 252
Left = 1560
TabIndex = 13
Top = 600
Width = 1332
End
End
Begin MSComDlg.CommonDialog CommonDialog1
Left = 240
Top = -120
_ExtentX = 847
_ExtentY = 847
_Version = 393216
CancelError = -1 'True
End
End
End
Attribute VB_Name = "FrmGrid"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Option Explicit
Private Declare Function PtInRegion Lib "gdi32" (ByVal hRgn As Long, ByVal X As Long, ByVal Y As Long) As Long
Private Declare Function CreatePolygonRgn Lib "gdi32" (lpPoint As POINTAPI, ByVal nCount As Long, ByVal nPolyFillMode As Long) As Long
Private Declare Function DeleteObject Lib "gdi32" (ByVal hObject As Long) As Long
Private Type POINTAPI
X As Long
Y As Long
End Type
Dim TheContouPath As String, TheContouFile As String, TheInstallPath As String
Dim Xmin0 As Double, Xmax0 As Double, Ymin0 As Double, Ymax0 As Double
Dim XminT As Double, XmaxT As Double, YminT As Double, YmaxT As Double
Dim DataType As Integer
Dim nDec As Integer, FMT As String
Dim Xmin As Double, Xmax As Double, DX As Double, NX As Integer
Dim Ymin As Double, Ymax As Double, DY As Double, NY As Integer
Dim Vmin As Single, Vmax As Single
Dim Xcontou() As Double, Ycontou() As Double, Zcontou() As Double, NContou As Long
Dim BorderX() As Single, BorderY() As Single, nBorder As Integer
Dim nSJX As Integer, ID1() As Integer, ID2() As Integer, ID3() As Integer
'按方位取点加权 weighting of measured points by azimuth *c
' NY : integer the number of lines in Y direction
' NX : integer the number of lines in X direction
' Ncontou : integer the number of original random data
' n0 : integer the number of azmuths in each 90 degree field
' Xmin : real X coordinater at the left-down of the lattice
' Ymin : real Y coordinater at the left-down of the lattice
' dx : real the interval of x
' dy : real the interval of y
' a : real (Ncontou,3) an array which store the original random data
' output formal parameters :
' s : real (NY,NX) an array which store the final data after latticized
'
Private Sub GridAZimuth(NY As Integer, NX As Integer, N0 As Integer, Xmin As Double, Ymin As Double, DX As Double, DY As Double, Xcontou() As Double, Ycontou() As Double, Zcontou() As Double, NContou As Long, Zgrid() As Double, bZgrid() As Boolean, Dr0 As Double)
Dim A1(0 To 9) As Double, B1(0 To 39) As Double, B2(0 To 39) As Double, Cn(0 To 39) As Double, Cp(0 To 39) As Double
Dim P As Integer, Pa As Double, C As Double, N2 As Integer, N3 As Integer, N4 As Integer
Dim Umi As Double, Vmi As Double, W As Double, X0 As Double, U As Double, V As Double, X As Double, Y As Double, Z As Double
Dim I As Integer, J As Integer, K As Integer, L As Integer
Dim I1 As Integer, J0 As Integer, J1 As Integer, L1 As Integer
Dim Y0 As Double, A() As Double
Dim XX0 As Double, YY0 As Double
C = 90# / N0 * 0.0174533
N4 = 4 * N0
N2 = 2 * N0
N3 = 3 * N0
Umi = Xmin / DX
Vmi = Ymin / DX
ReDim A(0 To NContou, 1 To 2)
For I = 0 To NContou
A(I, 1) = Xcontou(I) / DX
A(I, 2) = Ycontou(I) / DX
Next I
For I = 1 To N0 - 1
A1(I) = Tan(C * I)
Next I
A1(N0) = 90000000000#
For I = 1 To NY
Y0 = Vmi + (I - 1) * DY / DX
For J = 1 To NX
For J0 = 1 To 4 * N0
B1(J0) = 10000#
B2(J0) = 0#
Cp(J0) = 1#
Next J0
X0 = Umi + J - 1#
For K = 0 To NContou
X = A(K, 1) - X0
Y = A(K, 2) - Y0
Z = Zcontou(K)
If (X >= 0 And Y >= 0) Then
P = 0
ElseIf (X = 0 And Y >= 0) Then
P = N0
ElseIf (X = 0 And Y = 0) Then
P = N2
ElseIf (X >= 0 And Y = 0) Then
P = N3
End If
If (X = 0) Then
U = 1000000#
Else
If (N0 = 1) Then
U = 1#
Else
U = Abs(Y / X)
End If
End If
V = X * X + Y * Y
For L = 1 To N0
If (U A1(L)) Then
L1 = L + P
If (V B1(L1)) Then
B1(L1) = V
B2(L1) = Z
End If
Exit For
End If
Next L
Next K
Pa = 0#
For I1 = 1 To N4
If (B1(I1) > Pa) Then Pa = B1(I1)
Next I1
For I1 = 1 To N4
For J1 = 1 To N4
If (J1 <> I1) Then Cp(I1) = Cp(I1) * B1(J1) / Pa
Next J1
Next I1
W = 0#
For I1 = 1 To N4
W = W + Cp(I1)
Next I1
For I1 = 1 To N4
Cn(I1) = Cp(I1) / W
Next I1
Zgrid(J, I) = 0#
For I1 = 1 To N4
Zgrid(J, I) = Zgrid(J, I) + Cn(I1) * B2(I1)
Next I1
bZgrid(J, I) = True
YY0 = YminT + (I - 1) * DY
XX0 = XminT + (J - 1) * DX
PictureGrid.Circle (XX0, YY0), Dr0
DoEvents
Next J
Next I
End Sub
'趋势面拟合Fitting of trend surface
'input formal parameters
'NY : integer the number of lines in Y directioni
'NX : integer the number of points in X direction
'Ncontou : integer the number of original random data
'n0 : integer the number of times of trend surface
'xmin : real X coordinats at the left-down of the lattice
'ymin : real Y coordinats at the left-down of the lattice
'dx : real the interval of X
'dy : real the interval of Y
'd : real (Ncontou,3) an array which store the original random data
'output formal parameters
'd1 : real (Ncontou,3) an array which store residual erroes
'b : real (NY,NX) an array which store the trend surface of the lattice
Private Sub GridTrendSurface(NY As Integer, NX As Integer, N0 As Integer, Xmin As Double, Ymin As Double, DX As Double, DY As Double, Xcontou() As Double, Ycontou() As Double, Zcontou() As Double, NContou As Long, Zgrid() As Double, bZgrid() As Boolean, Dr0 As Double)
'D1(Ncontou, 3), B(NY, NX)
Dim H() As Double, S() As Double
Dim Nt As Integer, I As Integer, J As Integer, K As Integer, I0 As Integer, J0 As Integer
Dim X As Double, Y As Double, Z As Double, Sma As Double, B1 As Double
Dim E As Double, F As Double
Dim I9 As Integer, J9 As Integer
Dim XX0 As Double, YY0 As Double
Nt = 1
For I = 1 To N0
Nt = Nt + I + 1
Next I
ReDim S(1 To Nt, 1 To Nt + 1), H(1 To Nt)
For I = 1 To Nt
For J = 1 To Nt + 1
S(I, J) = 0#
Next J
Next I
For I = 0 To NContou
X = Xcontou(I)
Y = Ycontou(I)
Z = Zcontou(I)
X = ((X - Xmin) / DX) / (NX - 1)
Y = ((Y - Ymin) / DY) / (NY - 1)
For J = 1 To Nt
Call TSXY(X, Y, E, F, J)
H(J) = E * F
Next J
For I0 = 1 To Nt
For J0 = 1 To Nt
S(I0, J0) = H(I0) * H(J0) + S(I0, J0)
Next J0
S(I0, Nt + 1) = H(I0) * Z + S(I0, Nt + 1)
Next I0
Next I
For K = 1 To Nt
Sma = 0
For I = K To Nt
For J = 1 To Nt
If (Abs(Sma) Abs(S(I, J))) Then
Sma = S(I, J)
I9 = I
J9 = J
End If
Next J
Next I
'?????????????????if(sma.eq.0) goto 100
For J = 1 To Nt + 1
B1 = S(I9, J) / Sma
S(I9, J) = S(K, J)
S(K, J) = B1
Next J
For I = 1 To Nt
If (I <> K) Then
B1 = -S(I, J9)
For J = 1 To Nt + 1
S(I, J) = S(I, J) + B1 * S(K, J)
Next J
End If
Next I
Next K
For J = 1 To Nt
For I = 1 To Nt
If (S(I, J) >= 0.5) Then H(J) = S(I, Nt + 1)
Next I
Next J
For I = 1 To NY
Y = (I - 1) / (NY - 1)
For J = 1 To NX
X = (J - 1) / (NX - 1)
Zgrid(J, I) = 0#
For I0 = 1 To Nt
Call TSXY(X, Y, E, F, I0)
Zgrid(J, I) = H(I0) * E * F + Zgrid(J, I)
Next I0
bZgrid(J, I) = True
YY0 = YminT + (I - 1) * DY
XX0 = XminT + (J - 1) * DX
PictureGrid.Circle (XX0, YY0), Dr0
DoEvents
Next J
Next I
'For I = 1 To Ncontou
' X = ((D(I, 1) - Xmin) / dx) / (NX - 1)
' Y = ((D(I, 2) - Ymin) / dy) / (NY - 1)
' w = 0
' For I0 = 1 To Nt
' Call TSXY(X, Y, E, F, I0)
' w = H(I0) * E * F + w
' Next I0
' D1(I, 3) = D(I, 3) - w
' D1(I, 1) = D(I, 1)
' D1(I, 2) = D(I, 2)
'Next I
End Sub
'this subroutine is used to count coefficient of the linear equation
Private Sub TSXY(X As Double, Y As Double, E As Double, F As Double, I0 As Integer)
Dim K As Integer, L As Integer
K = 0
Do
If ((K + 1) * (K + 2) / 2 >= I0) Then Exit Do
K = K + 1
Loop
L = (K + 1) * (K + 2) / 2 - I0
If (X <> 0 Or L <> 0) Then
E = X ^ L
Else
E = 1
End If
If (Y <> 0 Or (K - L) <> 0) Then
F = Y ^ (K - L)
Else
F = 1
End If
End Sub
Private Sub CommandGridCacel_Click()
Unload Me
End Sub
Private Sub CommandGridOK_Click()
Dim iGrid As Integer
DX = Val(TextXStep.Text)
DY = Val(TextYStep.Text)
NX = Val(TextXNX.Text)
NY = Val(TextYNY.Text)
iGrid = ListGrid.ListIndex
Call DrawOld
If (iGrid = 3) Then
Call SJX(PictureGrid)
End If
Call Grid(iGrid, PictureGrid)
End Sub
Private Sub SJX(PictureSJX As PictureBox)
Dim R As Double, RT As Double
Dim J As Integer, K As Integer
Dim A1 As Double, B1 As Double, C1 As Double, Tt As Double
Dim AA As Double, BB As Double, CC As Double, CosC As Double
Dim bCheck As Boolean, bEQBD As Boolean
Dim Ni As Integer, Nj As Integer, Pi As Integer, Pj As Integer, Pk As Integer
Dim M1 As Integer, M2 As Integer
Dim X1 As Double, Y1 As Double, X2 As Double, Y2 As Double, X3 As Double, Y3 As Double
Dim bSJX As Boolean
'边界段数
Dim nBD As Integer
'用于判断边界环是否搜索过
Dim bBD() As Boolean
'用于判断数据点是否在边界环内
Dim bPoint() As Byte
'第i边界起点BD,第i边对应顶点
Dim BD() As Integer, BDij() As Integer
'第i边界的上一节点、下一节点
Dim Nlast As Integer, Nnext As Integer
ReDim ID1(0 To 2 * NContou), ID2(0 To 2 * NContou), ID3(0 To 2 * NContou)
'Begin生成第一个三角形
Pi = 0
'找出距第一点最近的点2
RT = 1E+20
For J = 1 To NContou
R = (Xcontou(J) - Xcontou(Pi)) ^ 2 + (Ycontou(J) - Ycontou(Pi)) ^ 2
If (R RT) Then
RT = R
Pj = J
End If
Next J
'找出第三点
Tt = 0
For J = 1 To NContou
If (J <> Pj) Then
AA = (Xcontou(Pj) - Xcontou(J)) ^ 2 + (Ycontou(Pj) - Ycontou(J)) ^ 2
BB = (Xcontou(Pi) - Xcontou(J)) ^ 2 + (Ycontou(Pi) - Ycontou(J)) ^ 2
CC = (Xcontou(Pi) - Xcontou(Pj)) ^ 2 + (Ycontou(Pi) - Ycontou(Pj)) ^ 2
CosC = 1# - (AA + BB - CC) / (2# * Sqr(AA * BB))
If (CosC > Tt + 0.00001) Then
Tt = CosC
Pk = J
End If
End If
Next J
nSJX = 0
ID1(nSJX) = Pi
ID2(nSJX) = Pj
ID3(nSJX) = Pk
'End生成第一个三角形
X1 = Xcontou(Pi)
X2 = Xcontou(Pj)
X3 = Xcontou(Pk)
Y1 = Ycontou(Pi)
Y2 = Ycontou(Pj)
Y3 = Ycontou(Pk)
PictureSJX.DrawMode = 13
PictureSJX.ForeColor = QBColor(12)
PictureSJX.Line (X1, Y1)-(X2, Y2)
PictureSJX.Line (X2, Y2)-(X3, Y3)
PictureSJX.Line (X3, Y3)-(X1, Y1)
DoEvents
'定义边界环数组
ReDim BD(0 To NContou), BDij(0 To NContou), bBD(0 To NContou), bPoint(0 To NContou)
For J = 0 To NContou
bBD(J) = False
bPoint(J) = 0
Next J
'生成三条边
nBD = 3
BD(1) = Pi
BDij(1) = Pk
bBD(1) = True
BD(2) = Pj
BDij(2) = Pi
bBD(2) = True
BD(3) = Pk
BDij(3) = Pj
bBD(3) = True
bPoint(Pi) = 1
bPoint(Pj) = 1
bPoint(Pk) = 1
Do
Ni = Ni + 1
bEQBD = False
For J = Ni To nBD
If (bBD(J) = True) Then
bEQBD = True
Exit For
End If
Next J
If (bEQBD = False) Then
For J = 1 To Ni - 1
If (bBD(J) = True) Then
bEQBD = True
Exit For
End If
Next J
End If
Ni = J
If (bEQBD = False) Then Exit Do
'第Ni边的终点节点
If (Ni = nBD) Then
Nj = 1
Else
Nj = Ni + 1
End If
'第Ni边的后邻
If (Ni = 1) Then
Nlast = nBD
Else
Nlast = Ni - 1
End If
'第Ni边的前邻
If (Nj = nBD) Then
Nnext = 1
Else
Nnext = Nj + 1
End If
'第Ni边的三个顶点
Pi = BD(Ni)
Pj = BD(Nj)
Pk = BDij(Ni)
'Begin找下一点
bCheck = False
If (Xcontou(Pj) = Xcontou(Pi)) Then
C1 = Xcontou(Pk) - Xcontou(Pi)
Else
A1 = (Ycontou(Pj) - Ycontou(Pi)) / (Xcontou(Pj) - Xcontou(Pi))
B1 = Ycontou(Pj) - A1 * Xcontou(Pj)
C1 = Ycontou(Pk) - Xcontou(Pk) * A1 - B1
End If
If (C1 = 0#) Then
M1 = -1
Else
M1 = 1
End If
Tt = 0#
For K = 1 To NContou
If (K = Pi Or K = Pj Or bPoint(K) = 2) Then
Else
If (Xcontou(Pj) = Xcontou(Pi)) Then
C1 = Xcontou(K) - Xcontou(Pi)
Else
C1 = Ycontou(K) - Xcontou(K) * A1 - B1
End If
If (C1 = 0#) Then
M2 = -1
Else
M2 = 1
End If
If (M1 <> M2) Then
AA = (Xcontou(Pj) - Xcontou(K)) ^ 2 + (Ycontou(Pj) - Ycontou(K)) ^ 2
BB = (Xcontou(Pi) - Xcontou(K)) ^ 2 + (Ycontou(Pi) - Ycontou(K)) ^ 2
CC = (Xcontou(Pi) - Xcontou(Pj)) ^ 2 + (Ycontou(Pi) - Ycontou(Pj)) ^ 2
CosC = 1# - (AA + BB - CC) / (2# * Sqr(AA * BB))
If (CosC > Tt + 0.00001) Then
Tt = CosC
Pk = K
bCheck = True
End If
End If
End If
Next K
'End找下一点
If (bCheck = True) Then '找到下一点
bSJX = True
'找前邻、后邻
If (Pk = BD(Nlast)) Then '后邻,删除Pi点
'修改后邻点参数
BDij(Nlast) = Pi
bBD(Nlast) = True
bPoint(Pk) = 1
For J = Ni + 1 To nBD
BD(J - 1) = BD(J)
BDij(J - 1) = BDij(J)
bBD(J - 1) = bBD(J)
Next J
nBD = nBD - 1 '减少一个边
bPoint(Pi) = 2 'Pi点位于圈内
ElseIf (Pk = BD(Nnext)) Then '前邻,删除Pj点
'修改Ni点参数
BDij(Ni) = Pj
bBD(Ni) = True
bPoint(Pk) = 1
For J = Nj + 1 To nBD
BD(J - 1) = BD(J)
BDij(J - 1) = BDij(J)
bBD(J - 1) = bBD(J)
Next J
nBD = nBD - 1 '减少一个边
bPoint(Pj) = 2 'Pj点位于圈内
ElseIf (bPoint(Pk) 1) Then '下一点不在边界环上,插入Pk点
'修改Ni点参数
BDij(Ni) = Pj
bBD(Ni) = True
bPoint(Pk) = 1
'数据后推
For J = nBD To Ni + 1 Step -1
BD(J + 1) = BD(J)
BDij(J + 1) = BDij(J)
bBD(J + 1) = bBD(J)
Next J
nBD = nBD + 1
'Pk点参数
BD(Ni + 1) = Pk
BDij(Ni + 1) = Pi
bBD(Ni + 1) = True
Ni = Ni + 1
Else
bSJX = False
End If
If (bSJX = True) Then
'生成新三角形
nSJX = nSJX + 1
ID1(nSJX) = Pi
ID2(nSJX) = Pj
ID3(nSJX) = Pk
X1 = Xcontou(Pi)
X2 = Xcontou(Pj)
X3 = Xcontou(Pk)
Y1 = Ycontou(Pi)
Y2 = Ycontou(Pj)
Y3 = Ycontou(Pk)
PictureSJX.Line (X1, Y1)-(X2, Y2)
PictureSJX.Line (X2, Y2)-(X3, Y3)
PictureSJX.Line (X3, Y3)-(X1, Y1)
DoEvents
End If
Else '没找到下一点,该边界不再搜索
bBD(Ni) = False
End If
Loop
'整理边界
nBorder = nBD - 1
ReDim BorderX(0 To nBorder + 1), BorderY(0 To nBorder + 1)
For J = 0 To nBorder
BorderX(J) = Xcontou(BD(J + 1))
BorderY(J) = Ycontou(BD(J + 1))
Next J
BorderX(nBorder + 1) = BorderX(0)
BorderY(nBorder + 1) = BorderY(0)
End Sub
Private Sub CommandContouMeshWang_Click()
Dim FalseTrue As Boolean
On Error Resume Next
CommonDialog1.DialogTitle = "等值线数据文件"
CommonDialog1.FileName = TheContouFile
CommonDialog1.Filter = "*.Txt;*.Dat|*.TXT;*.DAT"
CommonDialog1.InitDir = TheContouPath
CommonDialog1.FilterIndex = 0
CommonDialog1.ShowOpen
If (Err = 0) Then '打开文件
TheContouFile = CommonDialog1.DialogTitle
TheContouPath = CommonDialog1.FileName
Call ReadContouFile
If (DataType = 0) Then
FalseTrue = False
MsgBox "已经是网格化数据，无需网格化！", vbOKOnly, "关于网格化"
Else
FalseTrue = True
DX = Sqr((Xmax - Xmin) / (NContou + 1) * (Ymax - Ymin))
DY = DX
NX = Fix((Xmax - Xmin) / DX) + 1
NY = Fix((Ymax - Ymin) / DY) + 1
DX = (Xmax - Xmin) / (NX - 1)
DY = (Ymax - Ymin) / (NY - 1)
DX = Val(Format(DX, FMT))
DY = Val(Format(DY, FMT))
NX = Fix((Xmax - Xmin) / DX) + 1
NY = Fix((Ymax - Ymin) / DY) + 1
If ((NX - 1) * DX Xmax - Xmin) Then NX = NX + 1
If ((NY - 1) * DY Ymax - Ymin) Then NY = NY + 1
End If
LabelXmin.Enabled = FalseTrue
TextXmin.Enabled = FalseTrue
LabelXmax.Enabled = FalseTrue
TextXmax.Enabled = FalseTrue
LabelXStep.Enabled = FalseTrue
TextXStep.Enabled = FalseTrue
LabelXNX.Enabled = FalseTrue
TextXNX.Enabled = FalseTrue
LabelYmin.Enabled = FalseTrue
TextYmin.Enabled = FalseTrue
LabelYmax.Enabled = FalseTrue
TextYmax.Enabled = FalseTrue
LabelYStep.Enabled = FalseTrue
TextYStep.Enabled = FalseTrue
LabelYNY.Enabled = FalseTrue
TextYNY.Enabled = FalseTrue
LabelGrid.Enabled = FalseTrue
ListGrid.Enabled = FalseTrue
PictureGrid.Enabled = FalseTrue
CommandGridOK.Enabled = FalseTrue
TextXmin.Text = Xmin
TextXmax.Text = Xmax
TextYmin.Text = Ymin
TextYmax.Text = Ymax
TextXStep.Text = DX
TextYStep.Text = DY
TextXNX.Text = NX
TextYNY.Text = NY
ListGrid.ListIndex = 4
Call DrawOld
End If
End Sub
'数据网格化
Private Sub Grid(iGrid As Integer, PictureGrid As Object)
Dim CloseSJX(0 To 3) As POINTAPI, hRgnSJX As Long
Dim I As Integer, K As Integer, II As Integer
Dim IX As Integer, IY As Integer
Dim X As Double, Y As Double, Z As Double, bPtInRegion As Long
Dim X0 As Double, Y0 As Double, X1 As Double, Y1 As Double, Z1 As Double, X2 As Double, Y2 As Double, Z2 As Double, X3 As Double, Y3 As Double, Z3 As Double, X4 As Double, Y4 As Double, Z4 As Double
Dim Ds As Double, Ds1 As Double, Ds2 As Double, Ds3 As Double, GoTo1 As Integer
Dim Dr0 As Double, Delta As Double
Dim L1(0 To 2) As Double, L2(0 To 2) As Double
Dim RR0 As Double, R As Double, R0 As Double, R1 As Double, R2 As Double, R3 As Double
Dim RR(1 To 50) As Double, VV(1 To 50) As Double, V As Double, C As Double
Dim Fact As Double, X0c As Long, Y0c As Long
Dim TheOutFile As String, N0 As Integer
Dim AA() As Double, BB() As Double, XX() As Double
Dim Xt() As Double, Yt() As Double, Zt() As Double, Wt() As Double, Nt As Integer, S As Double
Dim Zgrid() As Double, bZgrid() As Boolean, bGrid As Boolean
If (Xmax0 - Xmin0 > Ymax0 - Ymin0) Then
Fact = 500 / (Xmax0 - Xmin0)
Else
Fact = 500 / (Ymax0 - Ymin0)
End If
PictureGrid.ForeColor = QBColor(2)
Dr0 = 4
Dr0 = Dr0 * (Xmax0 - Xmin0) / (PictureGrid.Width - 2 * Dr0)
PictureGrid.DrawMode = 13
ReDim Zgrid(1 To NX, 1 To NY), bZgrid(1 To NX, 1 To NY)
For IX = 1 To NX
For IY = 1 To NY
bZgrid(IX, IY) = False
Next IY
Next IX
If (iGrid = 3) Then
'Begin对三角形区域内网格点插值
For I = 0 To nSJX
X1 = Xcontou(ID1(I))
Y1 = Ycontou(ID1(I))
Z1 = Zcontou(ID1(I))
X2 = Xcontou(ID2(I))
Y2 = Ycontou(ID2(I))
Z2 = Zcontou(ID2(I))
X3 = Xcontou(ID3(I))
Y3 = Ycontou(ID3(I))
Z3 = Zcontou(ID3(I))
'建立三角形区域
CloseSJX(0).X = (X1 - Xmin0) * Fact
CloseSJX(0).Y = (Y1 - Ymin0) * Fact
CloseSJX(1).X = (X2 - Xmin0) * Fact
CloseSJX(1).Y = (Y2 - Ymin0) * Fact
CloseSJX(2).X = (X3 - Xmin0) * Fact
CloseSJX(2).Y = (Y3 - Ymin0) * Fact
CloseSJX(3).X = CloseSJX(0).X
CloseSJX(3).Y = CloseSJX(0).Y
hRgnSJX = CreatePolygonRgn(CloseSJX(1), 3, 2)
Y0 = YminT - DY
For IY = 1 To NY
Y0 = Y0 + DY
X0 = XminT - DX
For IX = 1 To NX
X0 = X0 + DX
If (bZgrid(IX, IY) = False) Then '该网格点未网格化
X0c = (X0 - Xmin0) * Fact
Y0c = (Y0 - Ymin0) * Fact
bPtInRegion = PtInRegion(hRgnSJX, X0c, Y0c)
If (bPtInRegion > 0) Then '在区域内
Select Case iGrid
Case 0 '三角平面插值
L1(0) = Y2 - Y1
L1(1) = X1 - X2
L1(2) = Y1 * (X2 - X1) - X1 * (Y2 - Y1)
L2(0) = Y3 - Y0
L2(1) = X0 - X3
L2(2) = Y0 * (X3 - X0) - X0 * (Y3 - Y0)
Delta = L1(0) * L2(1) - L2(0) * L1(1)
X = (L1(1) * L2(2) - L2(1) * L1(2)) / Delta
Y = (-L1(0) * L2(2) + L2(0) * L1(2)) / Delta
Ds1 = (X - X1) ^ 2 + (Y - Y1) ^ 2
Ds = (X2 - X1) ^ 2 + (Y2 - Y1) ^ 2
Z = Z1 + Ds1 / Ds * (Z2 - Z1)
Ds1 = (X0 - X3) ^ 2 + (Y0 - Y3) ^ 2
Ds = (X - X3) ^ 2 + (Y - Y3) ^ 2
Z = Z3 + Ds1 / Ds * (Z3 - Z)
Case 1 '双二次Lagrange插值
Call Lagrange(X1, Y1, Z1, X2, Y2, Z2, X3, Y3, Z3, X0, Y0, Z, 2)
Case 2 '双三次Lagrange插值
Call Lagrange(X1, Y1, Z1, X2, Y2, Z2, X3, Y3, Z3, X0, Y0, Z, 3)
Case 3 '距离平方反比插值
R1 = (X1 - X0) ^ 2 + (Y1 - Y0) ^ 2
R2 = (X2 - X0) ^ 2 + (Y2 - Y0) ^ 2
R3 = (X3 - X0) ^ 2 + (Y3 - Y0) ^ 2
Z = (Z1 / R1 + Z2 / R2 + Z3 / R3) / (1# / R1 + 1# / R2 + 1# / R3)
End Select
Zgrid(IX, IY) = Z
bZgrid(IX, IY) = True
PictureGrid.Circle (X0, Y0), Dr0
DoEvents
Else
'Begin判断是否在边线或顶点上
Ds1 = Sqr((X0 - X1) ^ 2 + (Y0 - Y1) ^ 2)
Ds2 = Sqr((X0 - X2) ^ 2 + (Y0 - Y2) ^ 2)
Ds = Sqr((X2 - X1) ^ 2 + (Y2 - Y1) ^ 2)
If (Abs(Ds - Ds1 - Ds2) 0.00001) Then '网格点为第一、二点边线上
Zgrid(IX, IY) = Z1 + Ds1 / Ds * (Z2 - Z1)
bZgrid(IX, IY) = True
PictureGrid.Circle (X0, Y0), Dr0
DoEvents
GoTo GoTo1
End If
Ds1 = Sqr((X0 - X2) ^ 2 + (Y0 - Y2) ^ 2)
Ds2 = Sqr((X0 - X3) ^ 2 + (Y0 - Y3) ^ 2)
Ds = Sqr((X3 - X2) ^ 2 + (Y3 - Y2) ^ 2)
If (Abs(Ds - Ds1 - Ds2) 0.00001) Then '网格点为第二、三点边线上
Zgrid(IX, IY) = Z2 + Ds1 / Ds * (Z3 - Z2)
bZgrid(IX, IY) = True
PictureGrid.Circle (X0, Y0), Dr0
DoEvents
GoTo GoTo1
End If
Ds1 = Sqr((X0 - X1) ^ 2 + (Y0 - Y1) ^ 2)
Ds2 = Sqr((X0 - X3) ^ 2 + (Y0 - Y3) ^ 2)
Ds = Sqr((X3 - X1) ^ 2 + (Y3 - Y1) ^ 2)
If (Abs(Ds - Ds1 - Ds2) 0.00001) Then '网格点为第一、三点边线上
Zgrid(IX, IY) = Z1 + Ds1 / Ds * (Z3 - Z1)
bZgrid(IX, IY) = True
PictureGrid.Circle (X0, Y0), Dr0
DoEvents
GoTo GoTo1
End If
'End判断是否在边线或顶点上
End If
GoTo1:
End If
Next IX
Next IY
DeleteObject (hRgnSJX)
Next I
'End对三角形区域内网格点插值
'Begin对三角形区域外网格点插值
PictureGrid.ForeColor = QBColor(12)
If (DX > DY) Then
R0 = 1.5 * DX
Else
R0 = 1.5 * DY
End If
RR0 = R0 * R0
Y0 = YminT - DY
For IY = 1 To NY
Y0 = Y0 + DY
X0 = XminT - DX
For IX = 1 To NX
X0 = X0 + DX
If (bZgrid(IX, IY) = False) Then
R0 = RR0
Do
II = 0
For K = 0 To NContou
R = (Xcontou(K) - X0) ^ 2 + (Ycontou(K) - Y0) ^ 2
If (R = R0) Then
II = II + 1
VV(II) = Zcontou(K)
RR(II) = R
End If
Next K
If (II > 0) Then
C = 0
For K = 1 To II
C = C + 1# / RR(K)
Next K
V = 0#
For K = 1 To II
V = V + VV(K) / RR(K)
Next K
V = V / C
Exit Do
End If
R0 = 1.5 * R0
bGrid = False
Loop
PictureGrid.Circle (X0, Y0), Dr0
DoEvents
Zgrid(IX, IY) = V
bZgrid(IX, IY) = True
End If
Next IX
Next IY
ElseIf (iGrid = 4) Then '按方位取点加权插值
Call GridAZimuth(NY, NX, 4, XminT, YminT, DX, DY, Xcontou, Ycontou, Zcontou, NContou, Zgrid, bZgrid, Dr0)
ElseIf (iGrid = 5) Then '趋势面拟合
N0 = Sqr(Sqr(NX * NY))
If (N0 4) Then N0 = 4
Call GridTrendSurface(NY, NX, N0, XminT, YminT, DX, DY, Xcontou, Ycontou, Zcontou, NContou, Zgrid, bZgrid, Dr0)
ElseIf (iGrid = 6) Then '加权最小二乘曲面拟合
ReDim BB(1 To 6), XX(1 To 6)
ReDim Xt(1 To 100), Yt(1 To 100), Zt(1 To 100), Wt(1 To 100)
PictureGrid.ForeColor = QBColor(12)
If (DX > DY) Then
R0 = 1.5 * DX
Else
R0 = 1.5 * DY
End If
RR0 = R0 * R0
Y0 = YminT - DY
For IY = 1 To NY
Y0 = Y0 + DY
X0 = XminT - DX
For IX = 1 To NX
X0 = X0 + DX
R0 = RR0
Do
Nt = 0
For K = 0 To NContou
R = (Xcontou(K) - X0) ^ 2 + (Ycontou(K) - Y0) ^ 2
If (R = 0.00001) Then
R = 0#
V = Zcontou(K)
Exit For
ElseIf (R = R0) Then
Nt = Nt + 1
Xt(Nt) = Xcontou(K)
Yt(Nt) = Ycontou(K)
Zt(Nt) = Zcontou(K)
Wt(Nt) = 1# / R
End If
Next K
If (R = 0.00001) Then
Exit Do
ElseIf (Nt > 6) Then
ReDim AA(1 To Nt, 1 To 6)
'Z
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Zt(K)
Next K
BB(1) = S
'ZX
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Zt(K) * Xt(K)
Next K
BB(2) = S
'ZY
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Zt(K) * Yt(K)
Next K
BB(3) = S
'ZXY
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Zt(K) * Xt(K) * Yt(K)
Next K
BB(4) = S
'ZXX
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Zt(K) * Xt(K) * Xt(K)
Next K
BB(5) = S
'ZYY
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Zt(K) * Yt(K) * Yt(K)
Next K
BB(6) = S
'1-----------
'AA(1,1)
S = 0#
For K = 1 To Nt
S = S + Wt(K)
Next K
AA(1, 1) = S
'AA(1,2)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K)
Next K
AA(1, 2) = S
'AA(1,3)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Yt(K)
Next K
AA(1, 3) = S
'AA(1,4)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Xt(K)
Next K
AA(1, 4) = S
'AA(1,5)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Yt(K)
Next K
AA(1, 5) = S
'AA(1,6)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Yt(K) * Yt(K)
Next K
AA(1, 6) = S
'2-------------
'AA(2,1)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K)
Next K
AA(2, 1) = S
'AA(2,2)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Xt(K)
Next K
AA(2, 2) = S
'AA(2,3)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Yt(K) * Xt(K)
Next K
AA(2, 3) = S
'AA(2,4)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * (Xt(K) ^ 3)
Next K
AA(2, 4) = S
'AA(2,5)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Xt(K) * Yt(K)
Next K
AA(2, 5) = S
'AA(2,6)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Yt(K) * Yt(K)
Next K
AA(2, 6) = S
'3----------------------
'AA(3,1)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Yt(K)
Next K
AA(3, 1) = S
'AA(3,2)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Yt(K)
Next K
AA(3, 2) = S
'AA(3,3)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Yt(K) * Yt(K)
Next K
AA(3, 3) = S
'AA(3,4)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Xt(K) * Yt(K)
Next K
AA(3, 4) = S
'AA(3,5)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Yt(K) * Yt(K)
Next K
AA(3, 5) = S
'AA(3,6)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * (Yt(K) ^ 3)
Next K
AA(3, 6) = S
'4----------------
'AA(4,1)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Yt(K) * Xt(K)
Next K
AA(4, 1) = S
'AA(4,2)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Xt(K) * Yt(K)
Next K
AA(4, 2) = S
'AA(4,3)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Yt(K) * Yt(K)
Next K
AA(4, 3) = S
'AA(4,4)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Yt(K) * (Xt(K) ^ 3)
Next K
AA(4, 4) = S
'AA(4,5)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * ((Xt(K) * Yt(K)) ^ 2)
Next K
AA(4, 5) = S
'AA(4,6)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * (Yt(K) ^ 3)
Next K
AA(4, 6) = S
'5-------------
'AA(5,1)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Xt(K)
Next K
AA(5, 1) = S
'AA(5,2)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * (Xt(K) ^ 3)
Next K
AA(5, 2) = S
'AA(5,3)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Xt(K) * Yt(K)
Next K
AA(5, 3) = S
'AA(5,4)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * (Xt(K) ^ 4)
Next K
AA(5, 4) = S
'AA(5,5)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Yt(K) * (Xt(K) ^ 3)
Next K
AA(5, 5) = S
'AA(5,6)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * ((Xt(K) * Yt(K)) ^ 2)
Next K
AA(5, 6) = S
'6--------------
'AA(6,1)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Yt(K) * Yt(K)
Next K
AA(6, 1) = S
'AA(6,2)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * Yt(K) * Yt(K)
Next K
AA(6, 2) = S
'AA(6,3)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * (Yt(K) ^ 3)
Next K
AA(6, 3) = S
'AA(6,4)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * ((Xt(K) * Yt(K)) ^ 2)
Next K
AA(6, 4) = S
'AA(6,5)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * Xt(K) * (Yt(K) ^ 3)
Next K
AA(6, 5) = S
'AA(6,6)
S = 0#
For K = 1 To Nt
S = S + Wt(K) * (Yt(K) ^ 4)
Next K
AA(6, 6) = S
Call Gauss(AA, BB, 6, XX, N0)
V = XX(1) + XX(2) * X0 + XX(3) * Y0 + XX(4) * X0 * X0 + XX(5) * X0 * Y0 + XX(6) * Y0 * Y0
Exit Do
End If
R0 = 1.5 * R0
bGrid = False
Loop
PictureGrid.Circle (X0, Y0), Dr0
DoEvents
Zgrid(IX, IY) = V
bZgrid(IX, IY) = True
Next IX
Next IY
End If
Vmin = 10000000000#
Vmax = -Vmin
For IX = 1 To NX
For IY = 1 To NY
If (Zgrid(IX, IY) Vmin) Then Vmin = Zgrid(IX, IY)
If (Zgrid(IX, IY) > Vmax) Then Vmax = Zgrid(IX, IY)
Next IY
Next IX
'End对三角形区域外网格点插值
I = InStr(TheContouPath, ".")
If (I > 0) Then
TheOutFile = Left(TheContouPath, I - 1) + Format(iGrid, "0") + ".GRD"
Else
TheOutFile = TheContouPath + Format(iGrid, "0") + ".GRD"
End If
Open TheOutFile For Output As #1
Print #1, "DSAA"
Print #1, Format(NX, "####0 "); Format(NY, "####0")
Print #1, Format(XminT, "###0.000 "); Format(XmaxT, "###0.000")
Print #1, Format(YminT, "###0.000 "); Format(YmaxT, "###0.000")
Print #1, Format(Vmin, "###0.0000 "); Format(Vmax, "###0.0000")
For IY = 1 To NY
For IX = 1 To NX - 1
Print #1, Format(Zgrid(IX, IY), "####0.0000 ");
Next IX
Print #1, Format(Zgrid(NX, IY), "####0.0")
Next IY
Close (1)
End Sub
'读绘等值线的数据
Private Sub ReadContouFile()
Dim IX As Integer, IY As Integer, II As Integer
Dim I As Integer, J As Integer, K As Integer
Dim I1 As Integer, L As Integer, Temp As String, DSAA_DSBB As String * 4, ValueTemp As Single
Dim N0 As Long, N As Long
Dim Xt As Double, Yt As Double
Dim Lat As Double, Lon As Double, Rou As Double
Dim VminTMP As Double, VmaxTMP As Double
Dim iModeOld As Integer, StrMax As String
Dim bSpace As Boolean
Dim lNX As Long, lNY As Long, DouValueTemp As Double, xStep As Double, yStep As Double, VminT As Double, VmaxT As Double
Open TheContouPath For Binary Access Read Lock Read As #1
Get #1, 1, DSAA_DSBB
Close (1)
If (DSAA_DSBB = "DSAA") Then
Open TheContouPath For Input As #1
Line Input #1, Temp
Input #1, NX, NY
Input #1, Xmin, Xmax
Input #1, Ymin, Ymax
Input #1, Vmin, Vmax
For IX = 1 To NX
Input #1, ValueTemp
Next IX
Line Input #1, Temp
Close (1)
If (Len(Trim(Temp)) 1) Then
bSpace = True
Else
bSpace = False
End If
Close (1)
DataType = 0
Open TheContouPath For Input As #1
Line Input #1, Temp
Input #1, NX, NY
Input #1, Xmin, Xmax
Input #1, Ymin, Ymax
Input #1, Vmin, Vmax
ReDim Zgrid(1 To NX, 1 To NY), bZgrid(1 To NX, 1 To NY)
N0 = NX * NY
ReDim Xcontou(0 To N0), Ycontou(0 To N0), Zcontou(0 To N0)
NContou = -1
DX = (Xmax - Xmin) / (NX - 1)
DY = (Ymax - Ymin) / (NY - 1)
Yt = Ymin - DY
If (bSpace = False) Then
For IY = 1 To NY
Yt = Yt + DY
Xt = Xmin - DX
For IX = 1 To NX
Xt = Xt + DX
NContou = NContou + 1
Input #1, Zcontou(NContou)
Xcontou(NContou) = Xt
Ycontou(NContou) = Yt
Zgrid(IX, IY) = Zcontou(NContou)
Next IX
Next IY
Else
For IY = 1 To NY
Yt = Yt + DY
Xt = Xmin - DX
For IX = 1 To NX
Xt = Xt + DX
NContou = NContou + 1
Input #1, Zcontou(NContou)
Xcontou(NContou) = Xt
Ycontou(NContou) = Yt
Zgrid(IX, IY) = Zcontou(NContou)
Next IX
Line Input #1, Temp
Next IY
End If
Close (1)
ElseIf (DSAA_DSBB = "DSBB") Then
DataType = 0
Open TheContouPath For Binary Access Read Lock Read As #1
Seek #1, 5
Get #1, , NX
Get #1, , NY
Get #1, , Xmin
Get #1, , Xmax
Get #1, , Ymin
Get #1, , Ymax
Get #1, , Vmin
Get #1, , Vmax
ReDim Zgrid(1 To NX, 1 To NY), bZgrid(1 To NX, 1 To NY)
N0 = NX
N0 = N0 * NY
ReDim Xcontou(0 To N0), Ycontou(0 To N0), Zcontou(0 To N0)
NContou = -1
DX = (Xmax - Xmin) / (NX - 1)
DY = (Ymax - Ymin) / (NY - 1)
Yt = Ymin - DY
For IY = 1 To NY
Yt = Yt + DY
Xt = Xmin - DX
For IX = 1 To NX
Xt = Xt + DX
NContou = NContou + 1
Get #1, , ValueTemp
Xcontou(NContou) = Xt
Ycontou(NContou) = Yt
Zcontou(NContou) = ValueTemp
Zgrid(IX, IY) = ValueTemp
Next IX
Next IY
Close (1)
ElseIf (DSAA_DSBB = "DSRB") Then
DataType = 0
Open TheContouPath For Binary Access Read Lock Read As #1
Seek #1, 17
Get #1, , lNY
Get #1, , lNY
Get #1, , lNX
Get #1, , XminT
Get #1, , YminT
Get #1, , xStep
Get #1, , yStep
Get #1, , VminT
Get #1, , VmaxT
NX = lNX
NY = lNY
DX = xStep
DY = yStep
Xmin = XminT
Xmax = Xmin + (NX - 1) * xStep
Ymin = YminT
Ymax = Ymin + (NY - 1) * yStep
Vmin = VminT
Vmax = VmaxT
ReDim Zgrid(1 To NX, 1 To NY), bZgrid(1 To NX, 1 To NY)
N0 = NX
N0 = N0 * NY
ReDim Xcontou(0 To N0), Ycontou(0 To N0), Zcontou(0 To N0)
Seek #1, 101
NContou = -1
Yt = Ymin - DY
For IY = 1 To NY
Yt = Yt + DY
Xt = Xmin - DX
For IX = 1 To NX
Xt = Xt + DX
NContou = NContou + 1
Get #1, , DouValueTemp
Xcontou(NContou) = Xt
Ycontou(NContou) = Yt
Zcontou(NContou) = DouValueTemp
Zgrid(IX, IY) = DouValueTemp
Next IX
Next IY
Close (1)
Else
'判断一行有几个数
Open TheContouPath For Input As #1
For I = 1 To 3
Line Input #1, Temp
Temp = Trim(Temp)
J = Len(Temp)
I1 = 2
K = 1
Do While I1 J
If (Mid(Temp, I1, 1) = " " Or Mid(Temp, I1, 1) = ",") Then
K = K + 1
For L = I1 + 1 To J
If (Mid(Temp, L, 1) = " " Or Mid(Temp, L, 1) = ",") Then
I1 = L
Else
I1 = I1 + 1
Exit For
End If
Next L
Else
I1 = I1 + 1
End If
Loop
If (K <> 3) Then Exit For
Next I
Close (1)
If (K <> 3) Then
MsgBox "非本程序识别格式!请按如下格式存放:" + Chr(10) + Chr(13) + "纬度,经度,观测值", vbOKOnly, "关于绘平面图等"
Exit Sub
End If
N0 = 1000
ReDim Xcontou(0 To N0), Ycontou(0 To N0), Zcontou(0 To N0)
N = -1
Open TheContouPath For Input As #1
Do While Not EOF(1)
N = N + 1
If (N > N0) Then
N0 = N0 + 100
ReDim Preserve Xcontou(0 To N0), Ycontou(0 To N0), Zcontou(0 To N0)
End If
Input #1, Ycontou(N), Xcontou(N), Zcontou(N)
Loop
NContou = N
Close (1)
Call SortXY(Xcontou, Ycontou, Zcontou, NContou)
'判断是否是网格数据
Call CheckGrid(Xcontou, Ycontou, NContou, NX, NY)
If (DataType = 0) Then
ReDim Zgrid(1 To NX, 1 To NY), bZgrid(1 To NX, 1 To NY)
II = -1
For IX = 1 To NX
For IY = 1 To NY
II = II + 1
Zgrid(IX, IY) = Zcontou(II)
Next IY
Next IX
End If
Xmin = Xcontou(0)
Xmax = Xcontou(0)
Ymin = Ycontou(0)
Ymax = Ycontou(0)
Vmin = Zcontou(0)
Vmax = Zcontou(0)
For I = 0 To NContou
If (Xcontou(I) Xmin) Then Xmin = Xcontou(I)
If (Xcontou(I) > Xmax) Then Xmax = Xcontou(I)
If (Ycontou(I) Ymin) Then Ymin = Ycontou(I)
If (Ycontou(I) > Ymax) Then Ymax = Ycontou(I)
If (Zcontou(I) Vmin) Then Vmin = Zcontou(I)
If (Zcontou(I) > Vmax) Then Vmax = Zcontou(I)
Next I
'Begin判断经纬度是否颠倒
If (Ymin >= -180 And Ymax = 180 And Xmin >= -90 And Xmax = 90) Then
ReDim Xcontou(0 To NContou), Ycontou(0 To NContou), Zcontou(0 To NContou)
Open TheContouPath For Input As #1
For I = 0 To NContou
Input #1, Xcontou(I), Ycontou(I), Zcontou(I)
Next I
Close (1)
Call SortXY(Xcontou, Ycontou, Zcontou, NContou)
'判断是否是网格数据
Call CheckGrid(Xcontou, Ycontou, NContou, NX, NY)
If (DataType = 0) Then
ReDim Zgrid(1 To NX, 1 To NY), bZgrid(1 To NX, 1 To NY)
II = -1
For IX = 1 To NX
For IY = 1 To NY
II = II + 1
Zgrid(IX, IY) = Zcontou(II)
Next IY
Next IX
End If
Xmin = Xcontou(0)
Xmax = Xcontou(0)
Ymin = Ycontou(0)
Ymax = Ycontou(0)
Vmin = Zcontou(0)
Vmax = Zcontou(0)
For I = 0 To NContou
If (Xcontou(I) Xmin) Then Xmin = Xcontou(I)
If (Xcontou(I) > Xmax) Then Xmax = Xcontou(I)
If (Ycontou(I) Ymin) Then Ymin = Ycontou(I)
If (Ycontou(I) > Ymax) Then Ymax = Ycontou(I)
If (Zcontou(I) Vmin) Then Vmin = Zcontou(I)
If (Zcontou(I) > Vmax) Then Vmax = Zcontou(I)
Next I
End If
'End判断经纬度是否颠倒
End If
'设置等值线参数
Call MaxnDec(Vmin, Vmax, nDec, StrMax, FMT)
Vmin = Val(Format(Vmin, FMT))
Vmax = Val(Format(Vmax, FMT))
Xmin0 = Xmin
Xmax0 = Xmax
Ymin0 = Ymin
Ymax0 = Ymax
End Sub
'判断是否为网格数据
Private Sub CheckGrid(X() As Double, Y() As Double, N As Long, NX As Integer, NY As Integer)
Dim IX As Long, IY As Long
Dim RR() As Double, VV() As Double, V As Double, C As Double
Dim Xtemp As Double, Ytemp As Double, Error1 As Integer
Dim I As Long
DataType = 0
'判断是否是网格数据
Xtemp = X(0)
Ytemp = Y(0)
NX = 0
NY = 0
DY = Y(2) - Y(1)
For I = 0 To N
If (X(I) <> Xtemp) Then Exit For
NY = NY + 1
Next I
If (NY 3) Then 'X方向数据小于3个,肯定非网格数据
DataType = 1
GoTo Error1
End If
NX = Fix((N + 1) / NY)
If (NY * NX = N + 1) Then '有可能是网格数据
DataType = 0
'判断X坐标是否等间距
DX = X(NY) - X(0)
For IX = 2 To NX - 1
Xtemp = X(IX * NY) - X((IX - 1) * NY)
If (Abs(Xtemp - DX) > 0.000001) Then
DataType = 1
GoTo Error1
End If
Next IX
'判断Y是否等间距
DY = Y(1) - Y(0)
For IY = 2 To NY - 1
Ytemp = Y(IY) - Y(IY - 1)
If (Abs(Ytemp - DY) > 0.000001) Then
DataType = 1
GoTo Error1
End If
Next IY
Else '肯定不是网格数据
DataType = 1
End If
Error1:
End Sub
'数据按X、Y排序
Private Sub SortXY(X() As Double, Y() As Double, Z() As Double, N As Long)
Dim N1 As Long, N2 As Long, X1 As Double
Dim I As Long, J As Long
'按X坐标排序
Call ShellSort(X, Y, Z, N, 0, N)
'按Y坐标排序
X1 = X(0)
N1 = 0
N2 = 0
For I = 1 To N
If (X(I) = X1) Then
N2 = N2 + 1
Else
If (N2 > N1) Then '相同X按Y坐标排序
Call ShellSort(Y, X, Z, N, N1, N2)
End If
X1 = X(I)
N1 = I
N2 = I
End If
Next I
If (N2 > N1) Then '相同X按Y坐标排序
Call ShellSort(Y, X, Z, N, N1, N2)
End If
'平均重合点
J = -1
X1 = Z(0)
N1 = 1
For I = 1 To N
If (Abs(X(I) - X(I - 1)) + Abs(Y(I) - Y(I - 1)) 0.00001) Then
X1 = X1 + Z(I)
N1 = N1 + 1
Else
J = J + 1
Z(J) = X1 / N1
X(J) = X(I - 1)
Y(J) = Y(I - 1)
X1 = Z(I)
N1 = 1
End If
Next I
J = J + 1
Z(J) = X1 / N1
X(J) = X(N)
Y(J) = Y(N)
N = J
End Sub
'排序子程序
Private Sub ShellSort(X() As Double, Y() As Double, Z() As Double, N As Long, N1 As Long, N2 As Long)
Dim B As Long, M As Long, L As Long, I As Long, J As Long, K As Long
Dim Xtemp As Double, Ytemp As Double, Ztemp As Double, jk As Long
B = N2 - N1 + 1
B = Log(B) / Log(2#)
M = B
L = 2 ^ M
For I = 1 To M
K = L - 1
L = L / 2
For J = K + N1 To N2
Xtemp = X(J)
Ytemp = Y(J)
Ztemp = Z(J)
jk = J - K
Do While (jk > N1 - 1 And X(jk) > Xtemp)
X(jk + K) = X(jk)
Y(jk + K) = Y(jk)
Z(jk + K) = Z(jk)
jk = jk - K
If (jk 0) Then Exit Do
Loop
X(jk + K) = Xtemp
Y(jk + K) = Ytemp
Z(jk + K) = Ztemp
Next J
Next I
End Sub
Private Sub MaxnDec(Ymin As Single, Ymax As Single, nDec As Integer, StrMax As String, FormatTMP As String)
Dim ValueMax As Single, Delta As Single, StrYmin As String, StrYmax As String
If (Abs(Ymax) > Abs(Ymin)) Then
ValueMax = Abs(Ymax)
Else
ValueMax = Abs(Ymin)
End If
Delta = Abs(Ymax - Ymin)
If (Delta 0.001 Or ValueMax 0.001) Then
nDec = 5
FormatTMP = "######0.00000"
StrYmin = Format(Ymin, "######0.00000")
StrYmax = Format(Ymax, "######0.00000")
ElseIf (Delta 0.01 Or ValueMax 0.01) Then
nDec = 4
FormatTMP = "######0.0000"
StrYmin = Format(Ymin, "######0.0000")
StrYmax = Format(Ymax, "######0.0000")
ElseIf (Delta 0.1 Or ValueMax 0.1) Then
nDec = 3
FormatTMP = "######0.000"
StrYmin = Format(Ymin, "######0.000")
StrYmax = Format(Ymax, "######0.000")
ElseIf (Delta 1# Or ValueMax 1#) Then
nDec = 2
FormatTMP = "######0.00"
StrYmin = Format(Ymin, "######0.00")
StrYmax = Format(Ymax, "######0.00")
ElseIf (Delta 10# Or ValueMax 10#) Then
nDec = 1
FormatTMP = "######0.0"
StrYmin = Format(Ymin, "######0.0")
StrYmax = Format(Ymax, "######0.0")
Else
FormatTMP = "######0"
StrYmin = Format(Ymin, "######0")
StrYmax = Format(Ymax, "######0")
nDec = 0
End If
StrMax = StrYmax
If (Len(StrYmin) > Len(StrMax)) Then StrMax = StrYmin
End Sub
'双二次Lagrange、双三次Lagrange插值网格化
Private Sub Lagrange(X1 As Double, Y1 As Double, Z1 As Double, X2 As Double, Y2 As Double, Z2 As Double, X3 As Double, Y3 As Double, Z3 As Double, XX As Double, YY As Double, ZZ As Double, b23 As Integer)
Dim S As Double, S1 As Double, S2 As Double, S3 As Double, L As Double, A As Double, B As Double, C As Double
Dim U1 As Double, U2 As Double, U3 As Double
Dim Q(1 To 10) As Double, Z(1 To 10) As Double
Dim I As Integer
'Begin求大三角形V1V2V3的面积S
A = Sqr((X2 - X1) ^ 2 + (Y2 - Y1) ^ 2)
B = Sqr((X3 - X2) ^ 2 + (Y3 - Y2) ^ 2)
C = Sqr((X3 - X1) ^ 2 + (Y3 - Y1) ^ 2)
L = (A + B + C) / 2
S = Sqr(L * (L - A) * (L - B) * (L - C))
'End求大三角形V1V2V3的面积S
'Begin求小三角形VV2V3、VV1V3、VV1V2的面积S1、S2、S3
A = Sqr((X2 - XX) ^ 2 + (Y2 - YY) ^ 2)
B = Sqr((X3 - X2) ^ 2 + (Y3 - Y2) ^ 2)
C = Sqr((X3 - XX) ^ 2 + (Y3 - YY) ^ 2)
L = (A + B + C) / 2
S1 = Sqr(L * (L - A) * (L - B) * (L - C))
A = Sqr((XX - X1) ^ 2 + (YY - Y1) ^ 2)
B = Sqr((X3 - XX) ^ 2 + (Y3 - YY) ^ 2)
C = Sqr((X3 - X1) ^ 2 + (Y3 - Y1) ^ 2)
L = (A + B + C) / 2
S2 = Sqr(L * (L - A) * (L - B) * (L - C))
A = Sqr((X2 - X1) ^ 2 + (Y2 - Y1) ^ 2)
B = Sqr((XX - X2) ^ 2 + (YY - Y2) ^ 2)
C = Sqr((XX - X1) ^ 2 + (YY - Y1) ^ 2)
L = (A + B + C) / 2
S3 = Sqr(L * (L - A) * (L - B) * (L - C))
'End求小三角形VV2V3、VV1V3、VV1V2的面积S1、S2、S3
'Begin求V相对于三角形V1V2V3的重心坐标u1、u2、u3
U1 = S1 / S
U2 = S2 / S
U3 = S3 / S
'End求V相对于三角形V1V2V3的重心坐标u1、u2、u3
If (b23 = 3) Then '三次
'Begin求双三次多项式函数类
Q(1) = 4.5 * U1 * (U1 - 1# / 3) * (U1 - 2# / 3)
Q(2) = 4.5 * U2 * (U2 - 1# / 3) * (U2 - 2# / 3)
Q(3) = 4.5 * U3 * (U3 - 1# / 3) * (U3 - 2# / 3)
Q(4) = 13.5 * U1 * U2 * (U1 - 1# / 3)
Q(5) = 13.5 * U1 * U2 * (U2 - 1# / 3)
Q(6) = 13.5 * U1 * U3 * (U1 - 1# / 3)
Q(7) = 13.5 * U2 * U3 * (U2 - 1# / 3)
Q(8) = 13.5 * U1 * U3 * (U3 - 1# / 3)
Q(9) = 13.5 * U2 * U3 * (U3 - 1# / 3)
Q(10) = 27# * U1 * U2 * U3
'End求双三次多项式函数类
'Begin求10个插值点z坐标
Z(1) = Z1
Z(2) = Z2
Z(3) = Z3
Z(4) = (Z(1) + 0.5 * Z(2)) / 1.5
Z(5) = (Z(1) + 2# * Z(2)) / 3#
Z(6) = (Z(1) + 0.5 * Z(3)) / 1.5
Z(7) = (Z(2) + 0.5 * Z(3)) / 1.5
Z(8) = (Z(1) + 2# * Z(3)) / 3#
Z(9) = (Z(2) + 2# * Z(3)) / 3#
Z(10) = (Z(4) + Z(5) + Z(6) + Z(7) + Z(8) + Z(9)) / 6
'End求10个插值点z坐标
'Begin求双三次Lagrange网格化值
ZZ = 0#
For I = 1 To 10
ZZ = ZZ + Z(I) * Q(I)
Next I
'End求双三次Lagrange网格化值
Else '二次
'Begin求双二次多项式函数类
Q(1) = U1 * (2 * U1 - 1)
Q(2) = U2 * (2 * U2 - 1)
Q(3) = U3 * (2 * U3 - 1)
Q(4) = 4# * U2 * U3
Q(5) = 4# * U1 * U3
Q(6) = 4# * U1 * U2
'End求双二次多项式函数类
'Begin求6个插值点z坐标
Z(1) = Z1
Z(2) = Z2
Z(3) = Z3
Z(4) = (Z(2) + Z(3)) / 2#
Z(5) = (Z(1) + Z(3)) / 2#
Z(6) = (Z(1) + Z(2)) / 2
'End求6个插值点z坐标
'Begin求双二次Lagrange网格化值
ZZ = 0#
For I = 1 To 6
ZZ = ZZ + Z(I) * Q(I)
Next I
'End求双二次Lagrange网格化值
End If
End Sub
Private Sub DrawOld()
Dim Dx0 As Double, Dy0 As Double, I As Integer, J As Integer
Dim X As Double, Y As Double, YmaxNew As Double, XmaxNew As Double
PictureGrid.Picture = LoadPicture()
PictureGrid.DrawWidth = 1
XminT = Val(TextXmin.Text)
XmaxT = Val(TextXmax.Text)
YminT = Val(TextYmin.Text)
YmaxT = Val(TextYmax.Text)
DX = Val(TextXStep.Text)
DY = Val(TextYStep.Text)
NX = Val(TextXNX.Text)
NY = Val(TextYNY.Text)
If (XmaxT XminT + DX * (NX - 1)) Then XmaxT = XminT + DX * (NX - 1)
If (YmaxT YminT + DY * (NY - 1)) Then YmaxT = YminT + DY * (NY - 1)
If (Xmin XminT) Then
Xmin0 = Xmin
Else
Xmin0 = XminT
End If
If (Ymin YminT) Then
Ymin0 = Ymin
Else
Ymin0 = YminT
End If
If (Xmax > XmaxT) Then
Xmax0 = Xmax
Else
Xmax0 = XmaxT
End If
If (Ymax > YmaxT) Then
Ymax0 = Ymax
Else
Ymax0 = YmaxT
End If
Dx0 = 30
Dy0 = 30
Dx0 = Dx0 * (Xmax0 - Xmin0) / (PictureGrid.Width - 2 * Dx0)
Dy0 = Dy0 * (Ymax0 - Ymin0) / (PictureGrid.Height - 2 * Dy0)
PictureGrid.ScaleLeft = Xmin0 - Dx0
PictureGrid.ScaleWidth = (Xmax0 + Dx0) - (Xmin0 - Dy0)
PictureGrid.ScaleTop = Ymax0 + Dy0
PictureGrid.ScaleHeight = (Ymin0 - Dy0) - (Ymax0 + Dy0)
PictureGrid.ForeColor = 0 ' QBColor(4)
Dx0 = 3
Dx0 = Dx0 * (XmaxT - XminT) / (PictureGrid.Width - 2 * Dx0)
PictureGrid.DrawMode = 13
For I = 0 To NContou
PictureGrid.Circle (Xcontou(I), Ycontou(I)), Dx0
Next I
PictureGrid.ForeColor = QBColor(8)
PictureGrid.DrawMode = 10
XmaxNew = XminT + (NX - 1) * DX
YmaxNew = YminT + (NY - 1) * DY
X = XminT - DX
For I = 1 To NX
X = X + DX
PictureGrid.Line (X, YminT)-(X, YmaxNew)
Next I
Y = YminT - DY
For I = 1 To NY
Y = Y + DY
PictureGrid.Line (XminT, Y)-(XmaxNew, Y)
Next I
PictureGrid.ForeColor = QBColor(0)
End Sub
Private Sub Form_Load()
TheInstallPath = App.Path
If (Right(TheInstallPath, 1) <> "\") Then
TheInstallPath = App.Path + "\"
End If
TheContouPath = TheInstallPath + "等值线数据\"
TheContouFile = ""
End Sub
Private Sub TextXmax_KeyDown(KeyCode As Integer, Shift As Integer)
If (KeyCode = 13) Then
Call DrawOld
End If
End Sub
Private Sub TextXmin_KeyDown(KeyCode As Integer, Shift As Integer)
If (KeyCode = 13) Then
Call DrawOld
End If
End Sub
Private Sub TextXNX_KeyDown(KeyCode As Integer, Shift As Integer)
If (KeyCode = 13) Then
Call DrawOld
End If
End Sub
Private Sub TextXStep_KeyDown(KeyCode As Integer, Shift As Integer)
If (KeyCode = 13) Then
Call DrawOld
End If
End Sub
Private Sub TextYmax_KeyDown(KeyCode As Integer, Shift As Integer)
If (KeyCode = 13) Then
Call DrawOld
End If
End Sub
Private Sub TextYmin_KeyDown(KeyCode As Integer, Shift As Integer)
If (KeyCode = 13) Then
Call DrawOld
End If
End Sub
Private Sub TextYNY_KeyDown(KeyCode As Integer, Shift As Integer)
If (KeyCode = 13) Then
Call DrawOld
End If
End Sub
Private Sub TextYStep_KeyDown(KeyCode As Integer, Shift As Integer)
If (KeyCode = 13) Then
Call DrawOld
End If
End Sub
Private Sub Gauss(A() As Double, B() As Double, N As Integer, X() As Double, L As Integer)
Dim JS() As Integer
Dim I As Integer, J As Integer, K As Integer
Dim T As Double, D As Double, IIS As Integer
ReDim JS(1 To N)
L = 1
For K = 1 To N - 1
D = 0#
For I = K To N
For J = K To N
If (Abs(A(I, J)) > D) Then
D = Abs(A(I, J))
JS(K) = J
IIS = I
End If
Next J
Next I
If (D + 1# = 1#) Then
L = 0
Else
If (JS(K) <> K) Then
For I = 1 To N
T = A(I, K)
A(I, K) = A(I, JS(K))
A(I, JS(K)) = T
Next I
End If
If (IIS <> K) Then
For J = K To N
T = A(K, J)
A(K, J) = A(IIS, J)
A(IIS, J) = T
Next J
T = B(K)
B(K) = B(IIS)
B(IIS) = T
End If
End If
If (L = 0) Then Exit Sub
For J = K + 1 To N
A(K, J) = A(K, J) / A(K, K)
Next J
B(K) = B(K) / A(K, K)
For I = K + 1 To N
For J = K + 1 To N
A(I, J) = A(I, J) - A(I, K) * A(K, J)
Next J
B(I) = B(I) - A(I, K) * B(K)
Next I
Next K
If (Abs(A(N, N)) + 1# = 1#) Then
L = 0
Exit Sub
End If
X(N) = B(N) / A(N, N)
For I = N - 1 To 1 Step -1
T = 0#
For J = I + 1 To N
T = T + A(I, J) * X(J)
Next J
X(I) = B(I) - T
Next I
JS(N) = N
For K = N To 1 Step -1
If (JS(K) <> K) Then
T = X(K)
X(K) = X(JS(K))
X(JS(K)) = T
End If
Next K
End Sub