import math
    import numpy as np
    import matplotlib.pyplot as plt
    
    nl = 0.030      # Manning の粗度係数（低水路）
    nh = 0.035      # Manning の粗度係数（高水敷）
    Q = 1500        # 流量 [cu.m/s]
    hd = 2.5        # 下流端の水深 [m]
    xs = (          # 下流端からの距離 [m]
           0,  500, 1000, 1200, 1800,
        2100, 2500, 3000, 3300, 3800)
    zs = (          # 河床高 [m]
        0.0, 0.5, 0.9, 0.8, 2.0,
        2.3, 3.0, 3.0, 3.5, 4.0)
    Bls = (         # 低水路幅 [m]
        300, 320, 280, 250, 300,
        300, 320, 350, 300, 250)
    dl = 1.5        # 低水路深 [m]
    Bh = 300        # 高水敷幅 [m]
    
    p5_3 = 5 / 3
    pa_3 = p5_3 + p5_3
    QQ = Q * Q
    B_nh = Bh / nh
    
    def compaund(hl):
        '''
        複断面水路の速度水頭と摩擦勾配を算定
        Args:
            hl (float): 低水路の水深 [m]
        Returns:
            float: 速度水頭 [m]
            float: 摩擦勾配
        Note:
            断面の特性量，流量および冪数はグローバル変数
        '''
        if hl > dl:
            hh = hl - dl
            Il = math.pow(hl, p5_3) * B_nl; Al = Bl * hl
            Ih = math.pow(hh, p5_3) * B_nh; Ah = Bh * hh
            Sf = QQ / (Il + Ih)**2
            avv_2g = (Il**2 / Al + Ih**2 / Ah) * Sf / \
                        (19.6 * (Al + Ah))
        else:
            Il = math.pow(hl, p5_3) * B_nl
            Sf = QQ / (Il * Il)
            avv_2g = QQ / (19.6 * (Bl * hl)**2)
        return avv_2g, Sf
    
    sz = len(xs)
    Hs  = np.empty(sz)  # 水位 [m]
    Qls = np.empty(sz)  # 低水路の流量 [cu.m/s]
    
    # 境界条件
    Bl = Bls[0]
    B_nl = Bl / nl
    hl = hd
    H = zs[0] + hd
    avv_2g, Sf = compaund(hd)
    E = H + avv_2g
    
    Hs[0]  = H
    Qls[0] = math.pow(hd, p5_3) * B_nl * math.sqrt(Sf)
    
    for i, x, z, Bl in zip(range(1, sz), xs[1:], zs[1:], Bls[1:]):
        B_nl = Bl / nl
        dx_2 = (x - xs[i-1]) / 2
        cnst = E + Sf * dx_2 - z
        while True:
            avv_2g, Sf = compaund(hl)
            Sfdx_2 = Sf * dx_2
            er = hl + avv_2g - Sfdx_2 - cnst
            if abs(er) < 1e-4: break
            hl -= er / (1 + (pa_3 * Sfdx_2 - 2 * avv_2g) / hl)
        H = hl + z 
        E = H + avv_2g
    
        Hs[i]  = H
        Qls[i] = math.pow(hl, p5_3) * B_nl * math.sqrt(Sf)
    
    plt.plot(xs, zs, label="$z_b$")
    plt.plot(xs, Hs, label="$H$")
    plt.xlabel("$x$ (m)")
    plt.ylabel("$z_b, H$ (m)")
    plt.legend()
    plt.grid()
    plt.show()
    
    print("|$x$<br>[m]|$z$<br>[m]|$H$<br>[m]|$Q_l$<br>[m<sup>3</sup>/s]|")
    print("|---|---|---|---|")
    for i in range(sz):
        print("|{:.0f}|{:.1f}|{:.3f}|{:.0f}|".format(
            xs[i], zs[i], Hs[i], Qls[i]
        ))