Attribute VB_Name = "modMethod"
'曲面_网格
Option Explicit
'近N点按距离加权平均法作曲面插值
'X：数据点X坐标数组
'Y：数据点Y坐标数组
'Z：数据点Z坐标(函数值)数组
'A：插值点X坐标
'B：插值点Y坐标
'F：插值点函数值
'N：取离插值点最近的数据点个数
Public Sub NDS(X() As Double, Y() As Double, Z() As Double, _
A As Double, B As Double, F As Double, N As Integer)
Dim I As Integer, J As Integer, IC As Integer, M As Integer
Dim S1 As Double, S2 As Double, D As Double
Dim DIS(1000) As Double
On Error GoTo errL
M = UBound(X, 1) '数据点数
'求各数据点与插值点的距离
For I = 1 To M
DIS(I) = (X(I) - A) ^ 2 + (Y(I) - B) ^ 2
Next I
S1 = 0#: S2 = 0#
For I = 1 To N
IC = 1
For J = 1 To M
If DIS(J) < DIS(IC) Then IC = J
Next J
If DIS(IC) < 0.0001 Then
'当距离很近时数据点函数值即为插值点函数值
F = Z(IC)
Exit Sub
End If
D = Sqr(DIS(IC))
S1 = S1 + Z(IC) / D
S2 = S2 + 1 / D
DIS(IC) = 10000000
Next I
F = S1 / S2
Exit Sub
errL:
MsgBox "不同的数据点有相同的X、Y坐标，造成除数为0"
End Sub
'高斯消去法解线性代数方程组
Public Sub GAU(N As Integer, A() As Double, X() As Double)
Dim I As Integer, J As Integer, K As Integer, C As Double
On Error Resume Next
For K = 1 To N - 1
For I = K + 1 To N
For J = K + 1 To N + 1
A(I, J) = A(I, J) - A(I, K) * A(K, J) / A(K, K)
Next J
Next I
Next K
X(N) = A(N, N + 1) / A(N, N)
For K = N - 1 To 1 Step -1
C = 0
For J = K + 1 To N
C = C + A(K, J) * X(J)
Next J
X(K) = (A(K, N + 1) - C) / A(K, K)
Next K
End Sub
'加权最小二乘法作曲面插值
'X：数据点X坐标数组
'Y：数据点Y坐标数组
'Z：数据点Z坐标(函数值)数组
'A：插值点X坐标
'B：插值点Y坐标
'F：插值点函数值
'K：加权函数类型，K=1，2，...，9
Public Sub WLSA(X() As Double, Y() As Double, Z() As Double, _
A As Double, B As Double, F As Double, K As Integer)
Dim I As Integer, J As Integer, M As Integer
Dim X1 As Double, X2 As Double, Y1 As Double, Y2 As Double
Dim TERM As Double, XT As Double, YT As Double
Dim XXT As Double, YYT As Double, XYT As Double, ZT As Double
Dim E(1 To 6, 1 To 7) As Double, U(1 To 6) As Double
M = UBound(X, 1) '数据点数
For I = 1 To M
X1 = X(I): Y1 = Y(I): X2 = X1 ^ 2: Y2 = Y1 ^ 2
Select Case K '确定加权方式
Case 0 'I型
TERM = 1 / ((X1 - A) ^ 2 + (Y1 - B) ^ 2 + 0.00001)
Case 1 'II型
TERM = 1 / ((X1 - A) ^ 2 + (Y1 - B) ^ 2 + 0.00001) ^ 4
Case 2 'III型
TERM = Exp((-(X1 - A) ^ 2 - (Y1 - B) ^ 2) * 0.1) / _
((X1 - A) ^ 2 + (Y1 - B) ^ 2 + 0.00001) '系数为0.1
Case 3 'III型
TERM = Exp((-(X1 - A) ^ 2 - (Y1 - B) ^ 2) * 0.01) / _
((X1 - A) ^ 2 + (Y1 - B) ^ 2 + 0.00001) '系数为0.01
Case 4 'III型
TERM = Exp((-(X1 - A) ^ 2 - (Y1 - B) ^ 2) * 0.001) / _
((X1 - A) ^ 2 + (Y1 - B) ^ 2 + 0.00001) '系数为0.001
Case 5 'III型
TERM = Exp((-(X1 - A) ^ 2 - (Y1 - B) ^ 2) * 0.0001) / _
((X1 - A) ^ 2 + (Y1 - B) ^ 2 + 0.00001) '系数为0.0001
Case 6 'III型
TERM = Exp((-(X1 - A) ^ 2 - (Y1 - B) ^ 2) * 0.00001) / _
((X1 - A) ^ 2 + (Y1 - B) ^ 2 + 0.00001) '系数为0.00001
Case 7 'III型
TERM = Exp((-(X1 - A) ^ 2 - (Y1 - B) ^ 2) * 0.000001) / _
((X1 - A) ^ 2 + (Y1 - B) ^ 2 + 0.00001) '系数为0.000001
End Select
XT = X1 * TERM: YT = Y1 * TERM
XXT = X2 * TERM: YYT = Y2 * TERM: XYT = X1 * YT
E(1, 1) = E(1, 1) + TERM: E(1, 2) = E(1, 2) + XT
E(1, 3) = E(1, 3) + YT: E(1, 4) = E(1, 4) + XYT
E(1, 5) = E(1, 5) + XXT: E(1, 6) = E(1, 6) + YYT
E(2, 4) = E(2, 4) + X2 * YT: E(2, 5) = E(2, 5) + X2 * XT
E(2, 6) = E(2, 6) + Y2 * XT: E(3, 6) = E(3, 6) + Y2 * YT
E(4, 4) = E(4, 4) + X2 * YYT: E(4, 5) = E(4, 5) + X2 * XYT
E(4, 6) = E(4, 6) + Y2 * XYT: E(5, 5) = E(5, 5) + X2 * XXT
E(6, 6) = E(6, 6) + Y2 * YYT
ZT = Z(I) * TERM
U(1) = U(1) + ZT: U(2) = U(2) + X1 * ZT
U(3) = U(3) + Y1 * ZT: U(4) = U(4) + X1 * Y1 * ZT
U(5) = U(5) + X2 * ZT: U(6) = U(6) + Y2 * ZT
Next I
E(2, 2) = E(1, 5): E(2, 3) = E(1, 4)
E(3, 3) = E(1, 6): E(3, 4) = E(2, 6)
E(3, 5) = E(2, 4): E(5, 6) = E(4, 4)
For I = 1 To 5
For J = I + 1 To 6
E(J, I) = E(I, J)
Next J
Next I
For I = 1 To 6
E(I, 7) = U(I)
Next I
GAU 6, E, U '用高斯消去法解线性代数方程组
'插值结果
F = U(1) + A * (U(2) + B * U(4) + A * U(5)) + B * (U(3) + B * U(6))
End Sub
'计算网格点的函数值（网格化时使用）
'X：数组，观测数据的X坐标
'Y：数组，观测数据的Y坐标
'B：数组，保存网格点的函数值
'LLL：方法
'NNN：近点按距离加权平均的点数（加权最小二乘拟合时无用）
'KKK：加权最小二乘拟合的类型（近点按距离加权平均时无用）
Public Sub GRID(X() As Double, Y() As Double, B() As Double, _
LLL As Integer, NNN As Integer, KKK As Integer)
Dim N1 As Integer, M As Integer, N As Integer
Dim I As Integer, J As Integer, I0 As Integer
Dim XX As Double, YY As Double, F As Double
Dim miX As Double, maX As Double, miY As Double, maY As Double
Dim DX As Double, DY As Double
N1 = UBound(X, 1) 'N1为观测点个数
M = UBound(B, 1) '网格的行数
N = UBound(B, 2) '网格的列数
miX = X(1): miY = Y(1): maX = X(1): maY = Y(1)
For I = 1 To N1
If X(I) < miX Then miX = X(I) '求观测值X坐标最小值
If Y(I) < miY Then miY = Y(I) '求观测值Y坐标最小值
If X(I) > maX Then maX = X(I) '求观测值X坐标最大值
If Y(I) > maY Then maY = Y(I) '求观测值Y坐标最大值
Next I
DX = (maX - miX) / (N - 1) '网格在X方向上的增量
DY = (maY - miY) / (M - 1) '网格在Y方向上的增量
For I = 1 To M
YY = miY + DY * (I - 1)
For J = 1 To N
XX = miX + DX * (J - 1)
If LLL = 0 Then Call NDS(X, Y, Z, XX, YY, F, NNN)
If LLL = 1 Then Call WLSA(X, Y, Z, XX, YY, F, KKK)
B(I, J) = F
Next J
Next I
End Sub
