リング状磁石の磁界分布

#概要
リング状磁石の磁界分布を計算することが目標です．
リング中心の磁界は，理論的に積分すれば簡単に求まるのですが，任意位置(x,y,z)に対する磁気力の強さはきれいな関数で求まりません．

```math
\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](https://ja.wolframalpha.com/input/?i=1%2Fsqrt%28a-b*cos%28x%29-c*sin%28x%29%29+dx%E3%81%AE%E7%A9%8D%E5%88%86)で計算した結果です．とてもめんどくさい関数になります．

今回は実際の磁石の寸法が決まっており，どんな感じになるかを知りたいので，数値シミュレーションで計算をしてみようと思います．

#計算方法
磁気力は

```math
f=k\frac{q1q2}{|\vec r|^3} \vec r
```
に従うとします．今回は，$f/q2$を各座標（x,y,z)で計算します．

リング全体（半径r,中心座標(0,0,0)，x-y平面に存在）で磁価q2を持つとすると線密度磁化mは

```math
2 \pi r m = q2
```
になります．
$d\theta$による磁気力$\vec {df}$は

```math
\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の三つのグラフを計算します．

```math
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$で積分します．数値積分で，単純に長方形積分をします．

```python
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

```

こんな感じに関数化できました．

#シミュレーション結果
設定

```math
r=1\\
kq1=1\\
m=1\\
z=1
```

![fx.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/864bd9d2-80d4-18ee-2197-4abe6862ce53.png)
![fz.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/a3927e7b-79e4-7f97-2831-f3d0b8899001.png)
![fx-fz.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/0b27b8b3-d14a-13bb-66f5-9f67bac97d68.png)

直感的に考察しても，結果と一致します．
![fx.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/af43a334-acd9-a7a7-c899-787c65e61152.png)
一つ意外だったのは，fzはリング内ではほとんど変化がないことです．

![fz.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/12e37b39-6e63-4e1f-e89d-52f58659a300.png)

zの位置を変えてシミュレーションを回すと次のようになります．

```math
z=3
```

![fx-3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/dc3c6d6c-e71f-f1a2-ec91-c47e23d61aaf.png)
![fz-3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/fea72d37-d7e0-bc1c-bdd9-cf7ca1c3f933.png)
![fx-fz-3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/cf5d8761-4166-c074-a319-3da52b05e229.png)

性質はあまり変わりませんが，山の大きさが変わりました．

#シミュレーションコード
```python
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を設定して再度シミュレーションします．

```math
\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)
```
となります．

```math
h=0.5\\
r=1\\
kq1=1\\
m=1\\
z=3
```
![fx.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/a61c7208-5262-e726-882c-1fed962b7fc6.png)

![fz.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/e4ee3c42-8413-fe29-fdd8-a2b5e4e634a1.png)
![x-z.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/405484/e60ef076-86a3-1fe9-bba4-1b354d989df6.png)

形状はあまり変わりませんでしたが，高さは大きく小さくなりました．


```python

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")

```