simulation.py 文件源码

def f_electric(P):
    n = len(P)
    a = 1e-9*np.array([p.a for p in P])
    apa = a[:, None] + a[None, :]


    x = 1e-9*np.array([p.x.flatten() for p in P])

    R = x[:, None, :] - x[None, :, :]
    r = np.sqrt(np.sum(R**2, 2) + 1e-100)
    R0 = R / (r**3)[:, :, None]

    q = np.array([float(p.charge) for p in P])
    const = qq**2 / (4.*np.pi*eperm*rpermw)

    QQ = q[:, None] * q[None, :]
    F = const * QQ[:, :, None] * R0
    #F[np.diag_indices_from(r)] = 0.
    tooclose = r <= apa
    R0i = R / (np.maximum(a[:, None], a[None, :])**3)[:, :, None]
    F[tooclose] = (const * QQ[:, :, None] * R0i)[tooclose]

    f = np.sum(F, 1).T.reshape(3*n, 1)

    return f