VBA实现的曲面插值算法

在日常工作中，除了碰到二维插值外，还经常需要进行三维插值，也就是通过(x,y)内插出z，这就需要使用曲面插值算法，这里水文工具集给出一个VBA实现的曲面插值算法，并给出具体实现代码，从而方便直接调用。

VBA曲面插值算法函数名称为GetBilinearInterpolation，具体代码如下：

'================================
' VB中进行曲面插值
'
' 
'================================
Public Function GetBilinearInterpolation( _
       xAxis As Variant, _
       yAxis As Variant, _
       zSurface As Variant, _
       xcoord As Double, _
       ycoord As Double) As Double

'first find 4 neighbouring points
    nx = UBound(xAxis, 1)
    ny = UBound(yAxis, 1)
    Dim lx As Single    'index of x coordinate of adjacent grid point to left of P
    Dim ux As Single    'index of x coordinate of adjacent grid point to right of P

    GetNeigbourIndices xAxis, xcoord, lx, ux

    Dim ly As Single  'index of y coordinate of adjacent grid point below P
    Dim uy As Single  'index of y coordinate of adjacent grid point above P

    GetNeigbourIndices yAxis, ycoord, ly, uy

    fQ11 = zSurface(lx, ly)
    fQ21 = zSurface(ux, ly)
    fQ12 = zSurface(lx, uy)
    fQ22 = zSurface(ux, uy)

    'if point exactly found on a node do not interpolate
    If ((lx = ux) And (ly = uy)) Then
        GetBilinearInterpolation = fQ11
        Exit Function
    End If

    x = xcoord
    y = ycoord

    x1 = xAxis(lx, 1)
    x2 = xAxis(ux, 1)
    y1 = yAxis(ly, 1)
    y2 = yAxis(uy, 1)

    'if xcoord lies exactly on an xAxis node do linear interpolation
    If (lx = ux) Then
        GetBilinearInterpolation = fQ11 + (fQ12 - fQ11) * (y - y1) / (y2 - y1)
        Exit Function
    End If

    'if ycoord lies exactly on an xAxis node do linear interpolation
    If (ly = uy) Then
        GetBilinearInterpolation = fQ11 + (fQ22 - fQ11) * (x - x1) / (x2 - x1)
        Exit Function
    End If

    fxy = fQ11 * (x2 - x) * (y2 - y)
    fxy = fxy + fQ21 * (x - x1) * (y2 - y)
    fxy = fxy + fQ12 * (x2 - x) * (y - y1)
    fxy = fxy + fQ22 * (x - x1) * (y - y1)
    fxy = fxy / ((x2 - x1) * (y2 - y1))

    GetBilinearInterpolation = fxy
End Function

Public Sub GetNeigbourIndices( _
       inArr As Variant, _
     x As Double, _
       ByRef lowerX As Single, _
       ByRef upperX As Single)
       
    N = UBound(inArr, 1)
    If x <= inArr(1, 1) Then
        lowerX = 1
        upperX = 1
    ElseIf x >= inArr(N, 1) Then
        lowerX = N
        upperX = N
    Else
        For i = 2 To N
            If x < inArr(i, 1) Then
                lowerX = i - 1
                upperX = i
                Exit For
            ElseIf x = inArr(i, 1) Then
                lowerX = i
                upperX = i
                Exit For
            End If
        Next i
    End If
End Sub