import math
    import numpy as np
    import matplotlib.pyplot as plt
    
    n = 0.02     # Mainning の粗度係数
    B = 200      # 川幅 [m]
    Q = 2000     # 流量 [cu.m/s]
    So = 1 / 100 # 河床勾配
    L = 5000     # 全長 [m]
    hL = 1.4     # 上流端の水深 [m]
    
    # 汎用変数
    p = 10 / 3
    qq = (Q / B)**2
    qq_2g = qq / 19.6
    nnqq = n * n * qq
    
    def main(dx):
    
        K = nnqq * dx / 2
    
        sz = L // abs(dx) + 1
        xs = np.linspace( L, 0, sz)
        zs = np.linspace(50, 0, sz)
        hs = np.empty(sz)
    
        # 境界条件
        h = hL; hs[0] = h
        Sfdx_2 = K / math.pow(h, p)
        cnst = zs[0] + h + qq_2g / (h * h) + Sfdx_2
    
        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)
            hs[i] = h
            cnst = z + h + vv_2g + Sfdx_2
    
        return xs, zs, hs 
    
    xs1, zs1, hs1 = main(dx=-500)
    plt.plot(xs1, hs1, label="$\Delta x$=500m")
    
    xs2, zs2, hs2 = main(dx=-100)
    plt.plot(xs2, hs2, label="$\Delta x$=100m")
    
    plt.xlabel("$x$ (m)")
    plt.ylabel("$h$ (m)")
    plt.legend()
    plt.grid()
    plt.show()
    
    print(f"|$x$ [m]|$z$ [m]|$\Delta x$ = 500 [m]<br>$h$ [m]|$\Delta x$ = 100 [m]<br>$h$ [m]|")
    print("|---|---|---|---|")
    for i, x, z, h in zip(range(len(xs2)), xs2, zs2, hs2):
        if i % 5 == 0:
            j = i // 5
            print(f"|{x:.0f}|{z:.1f}|{hs1[j]:.3f}|{h:.3f}|")
        else:
            print(f"|{x:.0f}|{z:.1f}||{h:.3f}|")