#概要
リング状磁石の磁界分布を計算することが目標です．
リング中心の磁界は，理論的に積分すれば簡単に求まるのですが，任意位置(x,y,z)に対する磁気力の強さはきれいな関数で求まりません．
\int \frac{1}{\sqrt{ a - b cos(x) - c sin(x)}} dx
=-(2 sec(tan^(-1)(b/c) + x) sqrt(a - b cos(x) - c sin(x)) sqrt((-c sqrt(b^2/c^2 + 1) + b cos(x) + c sin(x))/(a - c sqrt(b^2/c^2 + 1))) sqrt((c (sqrt(b^2/c^2 + 1) + sin(x)) + b cos(x))/(a + c sqrt(b^2/c^2 + 1))) F_1(1/2;1/2, 1/2;3/2;(a - b cos(x) - c sin(x))/(a + sqrt(b^2/c^2 + 1) c), (a - b cos(x) - c sin(x))/(a - sqrt(b^2/c^2 + 1) c)))/(c sqrt(b^2/c^2 + 1))
WolframAlphaで計算した結果です．とてもめんどくさい関数になります．
今回は実際の磁石の寸法が決まっており，どんな感じになるかを知りたいので，数値シミュレーションで計算をしてみようと思います．
#計算方法
磁気力は
f=k\frac{q1q2}{|\vec r|^3} \vec r
に従うとします．今回は，$f/q2$を各座標（x,y,z)で計算します．
リング全体（半径r,中心座標(0,0,0)，x-y平面に存在）で磁価q2を持つとすると線密度磁化mは
2 \pi r m = q2
になります．
$d\theta$による磁気力$\vec {df}$は
\vec {df} = k\frac{q1*mrd\theta}{|\vec r|^3} \vec r\\
\vec r = (x-rcos(\theta),y-rsin(\theta),z)
になるはずです．
$\vec f $は$\vec {df}$は$d\theta$で積分すれば計算できるはずです（ポテンシャル法もありますが，微分したくないのでベクトルでやってます）．
#シミュレーション
結果はz軸に対し線対称になるはずです．よって，y座標は0とします（y=0）．x-zでの結果のみで十分となります．
結果はベクトルで（x,y,z方向があるので），位置zを固定し，位置xに対するfx/q2,fy/q2,fz/q2の三つのグラフを計算します．
df_x = k\frac{q1*mr d\theta}{(x^2+z^2+r^2-2xr cos(\theta))^{3/2}} (x-rcos(\theta))\\
df_z = k\frac{q1*mr d\theta}{(x^2+z^2+r^2-2xr cos(\theta))^{3/2}} z
を$d\theta$で積分します．数値積分で，単純に長方形積分をします．
def func_fx(x,z):
dtheta = 0.01
theta = 0
fx=0
#fz=0
#パラメータ
r=1
kq=1
m=1 #線密度
while True:
if theta >= math.pi:
break
fx = fx + kq*m*r*dtheta/((x**2+z**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(x-r*math.cos(theta))
#fz = fz + kq*m*r*dtheta/((x**2+z**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(z)
theta = theta + dtheta
# print(fx)
# print(fz)
return fx
def func_fz(x,z):
dtheta = 0.01
theta = 0
#fx=0
fz=0
#パラメータ
r=1
kq=1
m=1 #線密度
while True:
if theta >= math.pi:
break
#fx = fx + kq*m*r*dtheta/((x**2+z**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(x-r*math.cos(theta))
fz = fz + kq*m*r*dtheta/((x**2+z**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(z)
theta = theta + dtheta
# print(fx)
# print(fz)
return fz
こんな感じに関数化できました．
#シミュレーション結果
設定
r=1\\
kq1=1\\
m=1\\
z=1
直感的に考察しても，結果と一致します．
一つ意外だったのは，fzはリング内ではほとんど変化がないことです．
zの位置を変えてシミュレーションを回すと次のようになります．
z=3
性質はあまり変わりませんが，山の大きさが変わりました．
#シミュレーションコード
import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def func_fx(x,z):
dtheta = 0.01
theta = 0
fx=0
#fz=0
#パラメータ
r=1
kq=1
m=1 #線密度
while True:
if theta >= math.pi:
break
fx = fx + kq*m*r*dtheta/((x**2+z**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(x-r*math.cos(theta))
#fz = fz + kq*m*r*dtheta/((x**2+z**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(z)
theta = theta + dtheta
# print(fx)
# print(fz)
return fx
def func_fz(x,z):
dtheta = 0.01
theta = 0
#fx=0
fz=0
#パラメータ
r=1
kq=1
m=1 #線密度
while True:
if theta >= math.pi:
break
#fx = fx + kq*m*r*dtheta/((x**2+z**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(x-r*math.cos(theta))
fz = fz + kq*m*r*dtheta/((x**2+z**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(z)
theta = theta + dtheta
# print(fx)
# print(fz)
return fz
#グラフ化
x = np.arange(0, 5.0, 0.01)
fx=[]
fz=[]
z_fix = 3
for i in x:
hoge = func_fx(i,z_fix)
fx= np.append(fx,hoge)
for i in x:
hoge = func_fz(i,z_fix)
fz= np.append(fz,hoge)
plt.plot(x,fx)
plt.plot(x,fz)
#3次元グラフ
#fx
x = np.arange(-3.0, 3.0, 0.1)
y = np.arange(-3.0, 3.0, 0.1)
X, Y = np.meshgrid(x, y)
#z_fix = 1
fx= func_fx(X**2+Y**2, z_fix)
fig = plt.figure()
ax = Axes3D(fig)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("fx(x, y)")
ax.plot_wireframe(X, Y, fx)
plt.show()
#fz
#z_fix = 1
fz= func_fz(X**2+Y**2, z_fix)
fig2 = plt.figure()
ax2 = Axes3D(fig2)
ax2.set_xlabel("x")
ax2.set_ylabel("y")
ax2.set_zlabel("fz(x, y)")
ax2.plot_wireframe(X, Y, fz)
plt.show()
#追記：NS
よく考えたら，NSの2極があるので，そこも考慮しなければいけなかったですね．
磁石長さhを設定して再度シミュレーションします．
\vec {df} = k\frac{q1*mrd\theta}{|\vec r|^3} \vec r - k\frac{q1*mrd\theta}{|\vec r2|^3} \vec r2\\
\vec r = (x-rcos(\theta),y-rsin(\theta),z)
\vec r2 = (x-rcos(\theta),y-rsin(\theta),z+h)
となります．
h=0.5\\
r=1\\
kq1=1\\
m=1\\
z=3
形状はあまり変わりませんでしたが，高さは大きく小さくなりました．
import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def func_fx(x,z):
dtheta = 0.01
theta = 0
fx=0
#fz=0
#パラメータ
r=1
kq=1
m=1 #線密度
h=0.5
while True:
if theta >= math.pi:
break
fx = fx + kq*m*r*dtheta/((x**2+z**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(x-r*math.cos(theta))
fx = fx - kq*m*r*dtheta/((x**2+(z+h)**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(x-r*math.cos(theta))
#fz = fz + kq*m*r*dtheta/((x**2+z**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(z)
theta = theta + dtheta
# print(fx)
# print(fz)
return fx
def func_fz(x,z):
dtheta = 0.01
theta = 0
#fx=0
fz=0
#パラメータ
r=1
kq=1
m=1 #線密度
h=0.5
while True:
if theta >= math.pi:
break
#fx = fx + kq*m*r*dtheta/((x**2+z**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(x-r*math.cos(theta))
fz = fz + kq*m*r*dtheta/((x**2+z**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(z)
fz = fz - kq*m*r*dtheta/((x**2+(z+h)**2+r**2-2*x*r*math.cos(theta))**(3/2)) *(z+h)
theta = theta + dtheta
# print(fx)
# print(fz)
return fz
#グラフ化
x = np.arange(0, 5.0, 0.01)
fx=[]
fz=[]
z_fix = 3
for i in x:
hoge = func_fx(i,z_fix)
fx= np.append(fx,hoge)
for i in x:
hoge = func_fz(i,z_fix)
fz= np.append(fz,hoge)
figxz = plt.figure()
plt.plot(x,fx)
plt.plot(x,fz)
figxz.savefig("x-z.png")
#3次元グラフ
#fx
x = np.arange(-3.0, 3.0, 0.1)
y = np.arange(-3.0, 3.0, 0.1)
X, Y = np.meshgrid(x, y)
#z_fix = 1
fx= func_fx(X**2+Y**2, z_fix)
fig = plt.figure()
ax = Axes3D(fig)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("fx(x, y)")
ax.plot_wireframe(X, Y, fx)
plt.show()
fig.savefig("fx.png")
#fz
#z_fix = 1
fz= func_fz(X**2+Y**2, z_fix)
fig2 = plt.figure()
ax2 = Axes3D(fig2)
ax2.set_xlabel("x")
ax2.set_ylabel("y")
ax2.set_zlabel("fz(x, y)")
ax2.plot_wireframe(X, Y, fz)
plt.show()
fig2.savefig("fz.png")