VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   4860
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   6525
   LinkTopic       =   "Form1"
   ScaleHeight     =   4860
   ScaleWidth      =   6525
   StartUpPosition =   3  '窗口缺省
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   495
      Left            =   4080
      TabIndex        =   0
      Top             =   3600
      Width           =   1455
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False

'四點(diǎn)偏心法，上下邊界條件均為水位
'T為時(shí)間間隔，JP為每小時(shí)里的時(shí)層數(shù)，CN為糙率，M為計(jì)算時(shí)段的小時(shí)數(shù)，
'p為斷面數(shù)，p1為間距數(shù)，
Const T = 180, JP = 20, P1 = 20, CN = 0.037, M = 36, P = 21

Dim Z0(P), Z1(P), Z2(P), W0(P), W1(P), W2(P), W3(P), ZP(2)
Private Sub Command1_Click()
'四點(diǎn)偏心法為兩層格式，定義兩個(gè)時(shí)層的流量Q，水位Z，河寬B，截面面積A，邊界斷面每個(gè)小時(shí)節(jié)點(diǎn)的水位ZB，斷面間距X
'定義斷面形狀的參數(shù)，Z0,Z1,Z2,W0,W1,W2,W3分別為斷面中兩個(gè)梯形的高和寬，ZP(2)是為了首末斷面插值而設(shè)，r為四點(diǎn)偏心的系數(shù)
IP = 1: R = 0.65
Dim Q(2, P), Z(2, P), B(2, P), A(2, P), ZB(2, M), X(P1)
'LL,MM,PP,RR分別為四點(diǎn)偏心格式里面的
'Qn-1=Ln+Mn*Qn
'Zn=Pn+Rn*Qn
'里的Ln，Mn，Pn，Rn
Dim LL(P), MM(P), PP(P), RR(P), ZZ(P)
Open "in.txt" For Input As #1
Open "out.txt" For Output As #2
'首末斷面第m個(gè)小時(shí)的水位
For i = 1 To 2
    For j = 1 To M
        Input #1, ZB(i, j)
    Next j
Next i
'令初始流量為100，如果格式穩(wěn)定，誤差將在計(jì)算中消除
For i = 1 To P
    Q(1, i) = 50
Next i
'內(nèi)插每個(gè)斷面初始水位
For i = 1 To P
    Z(1, i) = ZB(1, 1) + (ZB(2, 1) - ZB(1, 1)) * (i - 1) / P1
Next i
'間距
For i = 1 To P1
    Input #1, X(i)
Next i
'斷面形狀參數(shù)
For i = 1 To P
    Input #1, Z0(i), Z1(i), Z2(i)
Next i
For i = 1 To P
    Input #1, W0(i), W1(i), W2(i), W3(i)
Next
Dim it As Integer
'計(jì)算各個(gè)斷面面積和寬度
For i = 1 To P
    it = i
    Call A_and_B(A(1, i), B(1, i), Z(1, i), it)
Next i
'Print #2, "斷面";
For i = 1 To P
        Q(2, i) = Q(1, i)
        Z(2, i) = Z(1, i)
        A(2, i) = A(1, i)
        B(2, i) = B(1, i)
Next i
For i = 1 To 21
    Print #2, Tab(10 * i); "#" & (i);
Next i
Print #2,
For n = 1 To JP * M - JP
    'If dd > 0.02 Then n = n - 1
    If n Mod JP = 0 Then Print n / JP, "hour"

    '內(nèi)插該時(shí)層的水位值，因?yàn)檫吔鐥l件只有每個(gè)小時(shí)的水位值，而時(shí)間間隔為60s
    For i = 1 To 2
        ZP(i) = ZB(i, IP) + (n - (IP - 1) * JP) / JP * (ZB(i, IP + 1) - ZB(i, IP))
    Next i
    If n Mod JP = 0 Then IP = IP + 1
    '計(jì)算邊界斷面面積和河寬
    Z(2, 1) = ZP(1)
    Z(2, P) = ZP(2)
    'Call A_and_B(A(2, 1), B(2, 1), Z(2, 1), 1)
    'Call A_and_B(A(2, P), B(2, P), Z(2, P), P)
    '開(kāi)始第一趟掃描
    RR(1) = 0
    PP(1) = ZP(1)
100    For i = 1 To P1
        BM = 0.5 * ((1 - R) * (B(1, i) + B(1, i + 1)) + R * (B(2, i) + B(2, i + 1)))
        AM = 0.5 * ((1 - R) * (A(1, i) + A(1, i + 1)) + R * (A(2, i) + A(2, i + 1)))
        QM = 0.5 * ((1 - R) * (Q(1, i) + Q(1, i + 1)) + R * (Q(2, i) + Q(2, i + 1)))
        ZM = 0.5 * ((1 - R) * (Z(1, i) + Z(1, i + 1)) + R * (Z(2, i) + Z(2, i + 1)))
        C1 = 2 * R * T / X(i) / BM
        E1 = Z(1, i) + Z(1, i + 1) + ((1 - R) / R) * C1 * (Q(1, i) - Q(1, i + 1))
        A2 = 2 * R * T / X(i) * ((QM / AM) ^ 2 * BM - 9.8 * AM)
        C9 = T / X(i) * QM / AM
        C2 = 1 - 4 * R * C9
        D2 = 1 + 4 * R * C9
        Call A_and_B(F8, B8, ZM, i)
        Call A_and_B(F9, B9, ZM, i + 1)
        E2 = (1 - R) / R * A2 * (Z(1, i + 1) - Z(1, i)) + (1 - 4 * (1 - R) * C9) * Q(1, i + 1)
        E2 = E2 + (1 + 4 * (1 - R) * C9) * Q(1, i) + 2 * T * (QM / AM) ^ 2 * (F9 - F8) / X(i)
        E2 = E2 - 2 * T * 9.8 * CN * CN * QM * QM * BM / AM / AM / (AM / BM) ^ (1 / 3)
        Y1 = RR(i) - C1
        Y2 = A2 * RR(i) + C2
        Y3 = E1 - PP(i)
        Y4 = E2 - A2 * PP(i)
        Y5 = Y2 + A2 * Y1
        LL(i + 1) = (A2 * Y3 + Y4) / Y5
        MM(i + 1) = -(D2 + A2 * C1) / Y5
        PP(i + 1) = (Y2 * Y3 - Y1 * Y4) / Y5
        RR(i + 1) = (Y1 * D2 - Y2 * C1) / Y5
    Next i
    '開(kāi)始第二趟掃描
    ZZ(P) = ZP(2)
    Q(2, P) = (ZZ(P) - PP(P)) / RR(P)
    For i = P - 1 To 1 Step -1
        Q(2, i) = LL(i + 1) + MM(i + 1) * Q(2, i + 1)
        ZZ(i) = PP(i) + RR(i) * Q(2, i)
    Next i
    dd = Abs(Z(2, 1) - ZZ(1))
    For i = 2 To P
        If dd < Abs(Z(2, i) - ZZ(i)) Then dd = Abs(Z(2, i) - ZZ(i))
    Next i
    For i = 1 To P
        Z(2, i) = ZZ(i)
        Call A_and_B(A(2, i), B(2, i), Z(2, i), i)
    Next i
    If dd > 0.02 Then GoTo 100
  
    If n Mod JP = 0 Then
        Print #2, "Z";
        For i = 1 To P Step 2
            
            Print #2, Tab(10 * i); Format(Z(2, i), "##.##");
            
       
        Next
        Print #2,
        Print #2, "Q";
        For i = 2 To P1 Step 2
            
            Print #2, Tab(10 * i); Format(Q(2, i), "######");
        Next
        Print #2,
            
        
    End If
    For i = 1 To P
        Q(1, i) = Q(2, i)
        Z(1, i) = Z(2, i)
        A(1, i) = A(2, i)
        B(1, i) = B(2, i)
    Next i
Next n
Close #1
Close #2

        
        
    
End Sub
Private Sub A_and_B(A, B, Z, it)


ZH = Z - Z0(it)
If ZH <= Z1(it) - Z0(it) Then
   B = W0(it) + (W1(it) - W0(it)) / (Z1(it) - Z0(it)) * ZH
   A = (W0(it) + B) / 2 * ZH
Else
   B = W2(it) + (W3(it) - W2(it)) / (Z2(it) - Z1(it)) * (ZH - Z1(it) + Z0(it))
   A = (W2(it) + B) / 2 * (ZH - Z1(it) + Z0(it)) + (W0(it) + W1(it)) / 2 * (Z1(it) - Z0(it))
 
End If

End Sub