import math
    import numpy as np
    from particle import Particle
    
    n = 0.025       # Manning の祖度係数
    Q = 4000        # 流量 [cu.m/s]
    h = 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)
    Bs = (          # 川幅 [m]
        300, 320, 280, 250, 300,
        300, 320, 350, 300, 250)
    d = 0.01        # 粒径 [cm]
    
    sz = 10
    Hs = np.empty(sz)
    tss = np.empty(sz)
    qss = np.empty(sz)
    cs  = np.empty(sz)
    
    # 汎用変数
    p10_3 = 10 / 3
    QQ = Q * Q
    QQ_2g = QQ / 19.6
    nnQQ = n * n * QQ
    
    # Particle オブジェクト
    p = Particle(d)
    sgd = p.sgd
    wf = p.wf
    
    # 境界条件
    H = zs[0] + h
    BB = Bs[0]**2
    E = H + QQ_2g / (BB * h**2)
    Sf = nnQQ / (BB * math.pow(h, p10_3))
    
    Hs[0]  = H
    tss[0] = 98000 * h * Sf / sgd
    qss[0] = p.Itakura(tss[0])
    
    for i, x, z, B in zip(range(1, sz), xs[1:], zs[1:], Bs[1:]):
        BB = B * B
        qq_2g = QQ_2g / BB
    
        dx_2 = (x - xs[i-1]) / 2
        K = nnQQ / BB * dx_2
        cnst = E + Sf * dx_2 - z
        # 水深の初期推定値として直下の計算点の水深を当てる
        while True:
            vv_2g = qq_2g / (h * h)         # 速度水頭
            Sfdx_2 = K / math.pow(h, p10_3) # 損失水頭
            er = h + vv_2g - Sfdx_2 - cnst # 残差
            if abs(er) < 1e-4: break
            # ニュートン法による補正
            h -= er / (1 + (p10_3 * Sfdx_2 - 2 * vv_2g) / h)
        H = h + z
        E = H + vv_2g
        Sf = Sfdx_2 / dx_2
    
        Hs[i]  = H
        tss[i] = 98000 * h * Sf / sgd
        qss[i] = p.Itakura(tss[i])
    
    cs[-1] = qss[-1] / wf
    for i, x, B in zip(range(sz - 2, -1, -1), xs[-2::-1], Bs[-2::-1]):
        Q_Bdx = 100 * Q / (B * (xs[i+1] - x)) # cu.cm/s/sq.cm
        cs[i] = (qss[i] + cs[i+1] * Q_Bdx) / (wf + Q_Bdx)
    
    print("|$x$<br>[m]|$H$<br>[m]|$\\tau_*$|$u_*$<br>[cm/s]|$q_{su}$<br>[cm<sup>3</sup>/s/cm<sup>2</sup>]|$c$|$q_{su}/w_f$|")
    print("|---|---|---|---|---|---|---|")
    for x, H, z, ts, qs, c in zip(xs, Hs, zs, tss, qss, cs):
        print("|{:.0f}|{:.2f}|{:.3f}|{:.3f}|{:.3f}|{:.3f}|{:.3f}|".format(
            x, H, ts, math.sqrt(ts * sgd), qs, c, qs / wf
        ))