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