磁場分布シミュレータv0.1をpythonで

（追記 5/4）下記のプログラムは不完全です。平面にかかる電流素片のみからの磁場を足し合わせています。コイルを巻く方向に360°積分するのをサボっているので、正確ではありません。ちゃんとしたプログラムは、またいつか作ります。。。。 （5/5）作りました。

実はアイディアがあるのです。誰もまねしないだろうから、言いますが、ベクトルパルスマグネットを作ろうとしています。なぜかというと、世の中にそんなものはないし、面白そうだからです。ベクトルマグネットはありますが、ベクトルパルスマグネットはないです。ググっても出てこない。

３月末に思いついてから、割と良い考えだなーと思って今まで来ました。

電源の材料調達のために、モノコイルのLCR回路のシミュレーターはそれらしいのができました。

スプリットペアとなると、磁場の値や、分布が気になるところ。しかもスプリットペアを２つ組み合わせるとなるとなおのことです。

というわけでまた、ネット上のサンプルプログラムを改変して、それっぽい磁場分布シミュレータを作ってみました。もちろんver 0.1です。

計算の中身は電流が作る磁場のベクトル場を、どんどん足し合わせているだけの、しょぼいプログラムです。それをベクトル場を表示するmatplotlib.pyplot.quiverという関数を使って表示しています。

定量性はないけど、なんかそれっぽい！

学生のN田くんがきて、「ラインで切って、磁場分布がプロットできるようになると良いですね・・・では、さよなら」と、コメントを残して去って行った。。。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

"""
Magnetic field distribution simulator v0.1

A. Ikeda of UEC Tokyo
2 May 2022

References
[1]https://krajit.github.io/sympy/vectorFields/vectorFields.html

"""




### basic parameters for calculations ###
xrang = 25
yrang = 25
div = 60
limit = 60
x, y = np.meshgrid(np.linspace(-xrang, xrang, div), np.linspace(-yrang, yrang, div))



### definition of B field generated by a wire ###
def wireat_x(a, b):
    return -(y - b)/((x - a)**2 + (y - b)**2)
def wireat_y(a, b):
    return (x - a)/((x - a)**2 + (y - b)**2)

u = x - x
v = y - y



### parameters for the 1st split pair ###
bore = 10
gap = 20
hgap = gap/2
borer = bore/2
layer = 10
wind = 5



### calculating B fields of the 1st split pair ###
for la in np.arange(layer):
    for wi in np.arange(wind):
        u += wireat_x(wi + hgap, -(borer + la)) -wireat_x(wi + hgap, (borer + la))
        v += wireat_y(wi + hgap, -(borer + la)) -wireat_y(wi + hgap, (borer + la))

for la in np.arange(layer):
    for wi in np.arange(wind):
        u += wireat_x(-wi - hgap, -(borer + la)) -wireat_x(-wi - hgap, (borer + la))
        v += wireat_y(-wi - hgap, -(borer + la)) -wireat_y(-wi - hgap, (borer + la))



### parameters for the 2nd split pair ###
bore2 = 5
borer2 = bore2/2
gap2 = 10
hgap2 = gap2/2
layer2 = 5
wind2 = 5
cf = 1



### calculating B fields of the 1st split pair ###
for la in np.arange(layer2):
    for wi in np.arange(wind2):
        u += (wireat_x(-(borer2 + la), wi + hgap2) -wireat_x(borer2 + la, wi + hgap2))/cf
        v += (wireat_y(-(borer2 + la), wi + hgap2) -wireat_y(borer2 + la, wi + hgap2))/cf

for la in np.arange(layer2):
    for wi in np.arange(wind2):
        u += (-wireat_x((borer2 + la), -wi - hgap2) + wireat_x(- (borer2 + la), -wi - hgap2))/cf
        v += (-wireat_y((borer2 + la), -wi - hgap2) + wireat_y(- (borer2 + la), -wi - hgap2))/cf
        

ul = np.where((u > limit) | (u < -limit), 0, u)
vl = np.where((v > limit) | (v < -limit), 0, v)

plt.figure(figsize = (12,12))
plt.quiver(x,y,ul,vl) 
plt.savefig('out.pdf')
plt.show()