import math
    
    class Particle:
    
        d      = None # 粒径 [cm]
        sgd    = None # s: 水中比重
        sqr_sg = None # g: 重力加速度 [cm/sq.s]
        as_rs  = None # 板倉の式中の α* / ρs
        wf     = None # 沈降速度 [cm/s]
        tsc    = None # 無次元掃流力
        inv_sqrPi = 1 / math.sqrt(math.pi)
    
        def __init__(self, d, s=1.65):
            ''' コンストラクタ
    
            沈降速度 wf は Rubey の式、無次元限界掃流力 tsc は 岩垣の式 で算定
    
            Args:
                d (float): 粒径 [cm]
                s (float): 水中比重。既定値 1.65
            '''
            sg = s * 980
            nu6_d = 0.01 * 6 / d # 0.01 : 動粘性係数 [sq.cm/s]
            self.d = d
            self.sqrt_sg = math.sqrt(sg)
            self.sgd = sg * d
            self.as_rs = 0.14 / (1 + s)
            self.wf = math.sqrt(2 * self.sgd / 3 + nu6_d**2) - nu6_d
            if d > 0.303:
                self.tsc = 80.9 * d / self.sgd
            elif d > 0.118:
                self.tsc = 134.6 * math.pow(d, 31/22) / self.sgd
            elif d > 0.0565:
                self.tsc = 55 * d / self.sgd
            elif d > 0.0065:
                self.tsc = 8.41 * math.pow(d, 11/32) / self.sgd
            else:
                self.tsc = 226 * d / self.sgd
    
        def PWRI(self, sq_us, n = 0.025):
            ''' 佐藤・吉川・芦田の式(土研式)
    
            Args:
                sq_us (float): 摩擦速度の自乗値 [sq.cm/sq.s]
                n (float): マニングの粗度係数。既定値 0.025
            Return:
                float: 単位幅・単位時間当たりの掃流砂量 [cu.cm/s/cm]
            '''
            ts = sq_us / self.sgd
            F = 1 / (1 + 8 * (self.tsc / ts)**4)
            if n >= 0.025:
                return math.pow(sq_us, 1.5) * F * 0.623
            else:
                f = 0.623 * math.pow(40 * n, -3.5)
                return math.pow(sq_us, 1.5) * F * f
    
        def MPM(self, tse):
            ''' Meyer-Peter・Müller の式
    
            Args:
                tse (float): 無次元有効掃流力
            Return:
                float: 単位幅・単位時間当たりの掃流砂量 [cu.cm/s/cm]
            '''
            if tse > self.tsc:
                return 8 * math.pow((tse - self.tsc) * self.d, 1.5) * self.sqrt_sg
            else:
                return 0
    
        def Itakura(self, ts, useApprox=False):
            ''' 板倉の式
    
            Args:
                ts (float): 無次元掃流力
                useApprox (bool): 近似式を使用.既定は使用しない
            Return:
                float: 単位面積・単位時間当たりの巻上げ量 [cu.cm/s/sq.cm]
            '''
            if useApprox:
                qsu = 0.008 * (self.as_rs * (14 * ts - 0.9) / math.sqrt(ts) - self.wf)
            else:
                Bs_ts = 0.143 / ts 
                a = Bs_ts - 2 # 2: 1 / ηo
                w = self.inv_sqrPi * math.exp(-a * a) / math.erfc(a) + 2
                omega = w / Bs_ts - 1
                qsu = 0.008 * (self.as_rs * omega * math.sqrt(self.sgd / ts) - self.wf) 
            return max(qsu, 0) 
    
        def cb_c(self, h, Sf, kappa=0.4):
            ''' 基準点濃度と鉛直平均濃度の比
    
            Args:
                h (float): 水深 [m]
                Sf (float): 摩擦勾配
                kappa (float): カルマン定数.既定値 0.4
            Return:
                float: 基準点濃度 cb と鉛直平均濃度 c の比
            '''
            us = 100 * math.sqrt(9.8 * h * Sf) # cm/s
            beta = 6 * self.wf / (kappa * us)
            return beta / (1 - math.exp(-beta))
download