【Python】クーロン力（２）
- atliqu
- 2025年3月23日
元ネタは、「電磁気学問題集 寺本吉輝 著」。 見た目的には全然面白くないけど、せっかくコードかいたから保存しておく。
正六角形の６つの頂点の５つに電荷を置いた時のそれぞれの電荷が受けるクーロン力。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 |
## Python 3.12.2, numpy 2.2.3, matplotlib 3.10.1 import numpy as np import matplotlib.pyplot as plt # 3D座標設定 def coordinate_3d(axes, range_x, range_y, range_z): axes.set_xlabel("x", fontsize=16) axes.set_ylabel("y", fontsize=16) axes.set_zlabel("z", fontsize=16) axes.set_xlim(range_x[0], range_x[1]) axes.set_ylim(range_y[0], range_y[1]) axes.set_zlim(range_z[0], range_z[1]) #3D グリッドの描画 def grid_3d(axes, range_x, range_y, range_z): x, y, z = np.mgrid[range_x[0]: range_x[1]: 11j, range_y[0]: range_y[1]: 11j, range_z[0]: range_z[1]: 11j] axes.scatter(x, y, z, s=0.2, color='red') #3Dベクトルの描画 def visual_vector_3d(axes, loc, vector, color='red'): axes.quiver(loc[0], loc[1], loc[2], vector[0], vector[1], vector[2], color = color, length = 1, arrow_length_ratio = 0.3) axes.view_init(elev=0, azim=90) #電荷の描画 def charge_3d(axes, pos, col): axes.scatter(pos[0], pos[1], pos[2], s=10, color=col) #クーロン力 def calculate_coulombs_force_vectorized(charges, positions): num_charges = len(charges) forces = np.zeros_like(positions) for i in range(num_charges): for j in range(num_charges): if i!=j: r_vec = positions[j]-positions[i] r_mag = np.linalg.norm(r_vec) if r_mag>1e-12: r_hat = r_vec/r_mag #方向を表す単位ベクトル force_mag = -k_e*charges[i]*charges[j]/(r_mag**2) forces[i] += force_mag*r_hat return forces #クーロン定数 1/4πε0 k_e = 8.9875517873681764e9 #電荷量 charges = np.array([-1e-5, -1e-5, -1e-5, -1e-5, -1e-5]) #電荷の位置 positions = np.array([[0.5, np.sqrt(3)/2.0, 0.0], [-0.5, np.sqrt(3)/2.0, 0.0], [-1.0, 0.0, 0.0], [-0.5, -np.sqrt(3)/2.0, 0.0], [0.5, -np.sqrt(3)/2.0, 0.0]]) #クーロン力の計算 net_forces = calculate_coulombs_force_vectorized(charges, positions) #クーロン力の3D表示 fig = plt.figure(figsize=(10, 10)) ax = fig.add_subplot(111, projection='3d') coordinate_3d(ax, [-5, 5], [-5, 5], [-5, 5]) for i in range(len(positions)): visual_vector_3d(ax, positions[i], net_forces[i], 'red') charge_3d(ax, positions[i], 'red') grid_3d(ax, [-5, 5], [-5, 5], [-5, 5]) plt.show() |
まあ、こうなるわなぁ．． 3次元で表示する意味もない。