import math
    import numpy as np
    import matplotlib.pyplot as plt
    
    n = 0.02      # Manning の祖度係数
    So = 1 / 1000 # 河床勾配
    B = 200       # 川幅 [m]
    Q = 2000      # 流量 [cu.m/s]
    h = 5         # 下流端の水深 [m]
    z = 0         # 下流端の河床高 [m]
    dx = 500      # 計算距離刻み [m]
    L = 5000      # 流路長 [m]
    
    sz = L // dx + 1
    xs = np.linspace(0, L, sz)
    zs = xs * So
    Hs = np.empty(sz)
    
    # 汎用変数
    p = 10 / 3
    qq = (Q / B)**2
    qq_2g = qq / 19.6
    K = n * n * qq * dx / 2
    
    # 境界条件
    H = z + h; Hs[0] = H
    
    cnst = H + qq_2g / (h * h) + K / math.pow(h, p)
    for i, z in zip(range(1, sz), zs[1:]):
        cnst -= z
        # 水深の初期推定値として直下の計算点の水深を当てる
        while True:
            vv_2g = qq_2g / (h * h)        # 速度水頭
            Sfdx_2 = K / math.pow(h, p)    # 損失水頭
            er = h + vv_2g - Sfdx_2 - cnst # 残差
            if abs(er) < 1e-4: break
            # ニュートン法による補正
            h -= er / (1 + (p * Sfdx_2 - 2 * vv_2g) / h)
        H = h + z; Hs[i] = H
        cnst = H + vv_2g + Sfdx_2
    
    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]|")
    print("|---|---|---|")
    for i in range(sz):
        print(f"|{xs[i]:.0f}|{zs[i]:.1f}|{Hs[i]:.3f}|")