import numpy as np
    import matplotlib.pyplot as plt
    from particle import Particle
    
    sz = 50
    ds = np.power(10, np.linspace(-2, 0, sz))
    ts1 = np.empty(sz) # 浮遊限界上限
    ts2 = np.empty(sz) # 浮遊限界下限
    tsc = np.empty(sz) # 移動限界
    
    for i, d in enumerate(ds):
        p = Particle(d)
        wf = p.wf
        sgd = p.sgd
        ts1[i] = (1.67 * wf)**2 / sgd
        ts2[i] = (1.08 * wf)**2 / sgd
        tsc[i] = p.tsc
    
    ticks      = [ 0.01,  0.02,  0.05,  0.1,  0.2,  0.5,  1,  2,  5]
    ticklabels = ['0.01','0.02','0.05','0.1','0.2','0.5','1','2','5']
    ax = plt.gca()
    ax.set_xscale('log')  # x 軸を log スケールで
    ax.set_yscale('log')  # y 軸を log スケールで
    ax.set_xticks(ticks)
    ax.set_yticks(ticks)
    ax.set_xticklabels(ticklabels)
    ax.set_yticklabels(ticklabels)
    plt.plot(ds, ts1, label='$u_*/w_f=1.67$')
    plt.plot(ds, ts2, label='$u_*/w_f=1.08$')
    plt.plot(ds, tsc, label='critical')
    plt.xlabel('$d$ [cm]')
    plt.ylabel('$\\tau_*$')
    plt.xlim(np.min(ds), np.max(ds))
    plt.ylim(0.01, 5)
    plt.legend()
    plt.grid(which="both")
    plt.show()