常用水文预报算法和计算程序VB版

 

Public Sub GETtnk(R() As Double, q() As Double, DT As Integer, N As Double, k As Double)

'根据入流和出流时间序列，用矩法算出本场洪水的纳须单位线参数 n 和 k

'R      净雨过程
'Q      R相应的出流过程
'DT     计算时段长

Dim MI1 As Double, MI2  As Double, MO1 As Double, MO2  As Double, NI2 As Double, NO2 As Double

Dim X1 As Double, X2 As Double, X3 As Double
Dim i As Integer, XM As Double

Dim NR As Integer, NQ As Integer

    NR = UBound(R())
    NQ = UBound(q)
   
    i = 0
    X1 = 0
    X2 = 0
    X3 = 0

    For i = 1 To NR
       
        X1 = X1 + R(i) * (2 * i - 1)
        X2 = X2 + R(i)
        X3 = X3 + R(i) * (2 * i - 1) ^ 2
        
    Next
    MI1 = X1 / X2 * DT / 2
    MI2 = X3 / X2 * DT * DT / 4
    X1 = 0
    X2 = 0
    X3 = 0
    For i = 1 To NQ - 1
        XM = (q(i) + q(i + 1)) / 2
        X1 = X1 + XM * (2 * i - 1)
        X2 = X2 + XM
        X3 = X3 + XM * (2 * i - 1) ^ 2

    Next
   
    MO1 = X1 / X2 * DT / 2
    MO2 = X3 / X2 * DT * DT / 4
    NI2 = MI2 - MI1 ^ 2
    NO2 = MO2 - MO1 ^ 2
    N = (MO1 - MI1) ^ 2 / (NO2 - NI2)
    k = (NO2 - NI2) / (MO1 - MI1)
    
End Sub

Public Sub CUNIT(N As Double, KK As Double, DT As Integer, UH() As Double, MAX As Integer, M As Integer)
'//时段单位线计算
   
   
    Dim IT As Double, TT As Double, ni As Double
    Dim JC As Double
    Dim X1 As Double, X2 As Double, UNIT As Double, WC As Double
   
   
    N = Int(N)
    UH(1) = 0
    IT = 2
    TT = DT
10  X1 = 0
    X2 = 0
   
    For i = 1 To N
        ni = N - 1
        JC = 1
        If ni = 0 Then GoTo 30
            For j = 1 To ni
               JC = JC * j
            Next
'        End If
30      X1 = ((TT + DT) / KK) ^ ni / JC + X1
        X2 = (TT / KK) ^ ni / JC + X2
       
    Next
    X1 = X1 * Exp(DT / KK)
    UH(IT) = Exp(-TT / KK) * (X1 - X2)
    If ((UH(IT) < UH(IT - 1) And UH(IT) < 0.0005)) Then
        M = IT - 1
        UNIT = 0
        For i = 1 To M
            UNIT = UNIT - UH(i)
        Next
        WC = UNIT - 1
        If Abs(WC) >= 0.0005 Then
            UH(M) = UH(M) - WC
            If (UH(M) < 0.0005) Then M = M - 1
        End If
       
    Else
        TT = TT + DT
        IT = IT + 1
        GoTo 10
    End If

End Sub 

'  模块名：InterpModule.bas
'  函数名：INLagrn
'  功能：  用拉格朗日插值公式进行一元全区间不等距插值
'  参数：  n     - Integer型变量，给定结点的点数
'          x     - Double型一维数组，长度为n，存放给定的n个结点的值x(i)，要求x(1)<x(2)<...<x(n)
'          y     - Double型一维数组，长度为n，存放给定的n个结点的函数值y(i)，y(i) = f(x(i)), i=1,2,...,n
'          t   - Double型变量，存放指定的插值点的值
'  返回值：Double型，指定的查指点t的函数近似值f(t)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''