円環電流の作る磁束と磁場をPythonで計算する。

# 円環電流の磁場
円筒座標系を考える。
半径rcの円環電流が、中心軸をz軸として、z=zcの位置に置かれているとします。
このとき位置(r, z)における

* flux(r, z, rc, zc, ic) :磁束
* br(r, z, rc, zc, ic) :r方向の磁場
* bz(r, z, rc, zc, ic) :z方向の磁場

は次の関数で与えられます。
ただし、ここでいう磁束とは半径rでzの位置にある円の内部を通過する磁束のことです。
r, z, rc, zcの単位はメートル(m)です。icの単位はアンペア(A)です。

```python:ring_current.py
import numpy as np
from scipy.special import *
from scipy import constants as sc

u0 = sc.mu_0 # 真空の透磁率
e0 = sc.epsilon_0 # 真空の誘電率
c0 = sc.c # 光速
pi = sc.pi

def flux(r, z, rc, zc, ic):
    if 0 == r:
        return 0.0
    k= 4*r*rc/((r+rc)**2+(z-zc)**2)
    fx = (2*pi*r)*(u0*ic/pi/k**0.5)*(rc/r)**0.5
    fx *= (1-k/2)*ellipk(k)-ellipe(k)
    return fx

def br(r, z, rc, zc, ic):
    if 0 == r:
        return 0.0    
    k= 4*r*rc/((r+rc)**2+(z-zc)**2)
    d = (r-rc)**2+(z-zc)**2
    br = (u0*ic/4/pi)*(z-zc)/k**0.5/(r*rc)**0.5
    br *= 2*(2-k)*rc*ellipe(k)/d-k*ellipk(k)/r
    return br

def bz(r, z, rc, zc, ic):
    k= 4*r*rc/((r+rc)**2+(z-zc)**2)
    d = (r-rc)**2+(z-zc)**2
    bz = (u0*ic/4/pi)*k**0.5/(r*rc)**0.5
    bz *= ellipk(k)-(r**2-rc**2+(z-zc)**2)*ellipe(k)/d
    return bz
```

式の導出は下の通りです。

## Lorentzゲージ表式

```math
\begin{align}

&\boldsymbol{E} = -\frac{\partial \boldsymbol{A_{L}}}{\partial t}-\operatorname{grad} \phi_L\\

&\boldsymbol{B}=\operatorname{rot}\boldsymbol{A_{L}}\\

&\Bigl(\Delta-\frac{1}{c^2}\frac{\partial^{2}}{\partial t^{2}}\Bigr)\boldsymbol{A_{L}}=-\mu_{0}\boldsymbol{i}\\

&\Bigl(\Delta-\frac{1}{c^2}\frac{\partial^{2}}{\partial t^{2}}\Bigr)\phi_L=-\frac{1}{\epsilon_{0}}\rho\\

&\operatorname{div}\boldsymbol{A_{L}}+\frac{1}{c^2}\frac{\phi_L}{t}=0
\end{align}
```

静磁場の計算式では次の２式を用います。

```math
\begin{align}

&\Delta\boldsymbol{A}=-\mu_{0}\boldsymbol{i}\\

&\boldsymbol{B}=\operatorname{rot}\boldsymbol{A}
\end{align}
```
この式を積分形式で書き直します。

## 積分形式

```math
\begin{align}
&\boldsymbol{A(r)}=\frac{\mu_0}{4 \pi}\int d^{3}x \frac{\boldsymbol i(x)}{|\boldsymbol{r-x}|} \\
&\boldsymbol{B(r)}=\frac{\mu_0}{4 \pi}\int d^{3}x \frac{\boldsymbol i(x) \times (\boldsymbol{r-x})}{|\boldsymbol{r-x}|^3}
\end{align}
```

更に線積分の形に書き直します。

```math
\begin{align}
&\boldsymbol{A(r)}=\frac{\mu_{0}I}{4 \pi}\int d\boldsymbol{l} \frac{1}{|\boldsymbol{r-x}|} \\
&\boldsymbol{B(r)}=\frac{\mu_{0}I}{4 \pi}\int  \frac{d\boldsymbol{l} \times (\boldsymbol{r-x})}{|\boldsymbol{r-x}|^3}
\end{align}
```
## 線積分の実行

円環電流の位置ベクトル$\boldsymbol{x}$、静磁場を計算する場所の位置ベクトルを$\boldsymbol{r}$とします。

```math
\boldsymbol{x}=\begin{pmatrix}
R \cos\theta\\
R \sin\theta\\
0
\end{pmatrix}
```

```math
\boldsymbol{r}=\begin{pmatrix}
r\\
0\\
z
\end{pmatrix}
```

このとき次のような関係式が成り立ちます。

```math
d \boldsymbol{l}=\begin{pmatrix}
-R \sin\theta\\
R \cos\theta\\
0
\end{pmatrix}d\theta
```

```math
d \boldsymbol{l}\times(\boldsymbol{r-x})=\begin{pmatrix}
R z \cos\theta\\
R z \sin\theta\\
R(R-r\cos\theta)
\end{pmatrix}d\theta
```

```math
\begin{align}
|\boldsymbol{r-x}|&=\sqrt{(r-R\cos\theta)^{2}+(0-R\sin\theta)^{2}+z^{2}}\\
&=\sqrt{r^{2}+R^{2}+z^{2}-2Rr\cos\theta}
\end{align}
```

これらの関係式を用いて線積分を実行します。

```math
A_y(\boldsymbol{r})=\frac{\mu_{0}I}{\pi \sqrt{k}}\sqrt{\frac{R}{r}}\Bigl((1-\frac{k}{2})K(k)-E(k)\Bigr)
```

```math
B_x(r)=\frac{\mu_{0}I}{4 \pi}\frac{z}{\sqrt{k}\sqrt{Rr}}\Bigl(\frac{2(2-k)R E(k)}{d^{2}}-\frac{k K(k)}{r}\Bigr)

```

```math
B_{z}(r)=\frac{\mu_{0}I}{4 \pi}\frac{\sqrt{k}}{\sqrt{Rr}}\Bigl(K(k)-\frac{(r^{2}-R^{2}+z^2)E(k)}{d^{2}}\Bigr)
```
ただし、

```math
k=\frac{4 r R}{(r+R)^{2}+z^{2}}, d^{2}=(r-R)^{2}+z^{2}
```

他の成分$A_x$、$A_z$、$B_y$はゼロ。

## 円柱座標表示への書き直し

ここでの系は軸対象なので、円筒座標系での表記に変更して、且つコイルの位置を$Z$とすると下の式になります。

```math
A_\theta(\boldsymbol{r})=\frac{\mu_{0}I}{\pi \sqrt{k}}\sqrt{\frac{R}{r}}\Bigl((1-\frac{k}{2})K(k)-E(k)\Bigr)

```

```math
B_r(\boldsymbol{r})=\frac{\mu_{0}I}{4 \pi}\frac{z-Z}{\sqrt{k}\sqrt{Rr}}\Bigl(\frac{2(2-k)R E(k)}{d^{2}}-\frac{k K(k)}{r}\Bigr)

```

```math
B_{z}(\boldsymbol{r})=\frac{\mu_{0}I}{4 \pi}\frac{\sqrt{k}}{\sqrt{Rr}}\Bigl(K(k)-\frac{(r^{2}-R^{2}+(z-Z)^2)E(k)}{d^{2}}\Bigr)
```
ただし、

```math
k=\frac{4 r R}{(r+R)^{2}+(z-Z)^{2}}, d^{2}=(r-R)^{2}+(z-Z)^{2}
```

## 磁束の式
```math
\Phi=\int\boldsymbol{B}\cdot d\boldsymbol{S}=\int\operatorname{rot}\boldsymbol{A}\cdot d\boldsymbol{S}=\int \boldsymbol{A}\cdot d\boldsymbol{l}
```
という関係が成り立つので、下の式になります。

```math
\Phi=2\pi r \times (線積分方向の\boldsymbol{A}の長さ)
```