Biot-Savart式

 

まずはBiot-Savartの式から

Q:=matrix([q[x],q[y],q[z]])://電流モーメント
A:=matrix([a[x],a[y],a[z]])://電流位置を原点とするセンサーの位置ベクトル
B:=Symbol::mu[0]/(4*PI)*linalg::crossProduct(Q,A)/norm(A,2)^3;

matrix([[-((a[y]*q[z] - a[z]*q[y])*`μ`[0])/(4*PI*(abs(a[x])^2 + abs(a[y])^2 + abs(a[z])^2)^(3/2))], [((a[x]*q[z] - a[z]*q[x])*`μ`[0])/(4*PI*(abs(a[x])^2 + abs(a[y])^2 + abs(a[z])^2)^(3/2))], [-((a[x]*q[y] - a[y]*q[x])*`μ`[0])/(4*PI*(abs(a[x])^2 + abs(a[y])^2 + abs(a[z])^2)^(3/2))]])

電流の位置を(10,0,0)cm、大きさを(0,0,20)nA、センサーの位置をR(r[x],r[y],r[z])とします。

rule1:=[Symbol::mu[0]=4*PI*1E-7]:
rule2:=[a[x]=r[x]-10E-2,a[y]=r[y],a[z]=r[z]]:
rule3:=[q[x]=0,q[y]=0,q[z]=20E-9]:
E1:=subs(B,[op(rule1),op(rule2),op(rule3)]);

matrix([[-(2.0*10^(-15)*r[y])/(abs(r[x] - 0.1)^2 + abs(r[y])^2 + abs(r[z])^2)^(3/2)], [(0.0000001*(0.00000002*r[x] - 0.000000002))/(abs(r[x] - 0.1)^2 + abs(r[y])^2 + abs(r[z])^2)^(3/2)], [0]])

x軸、y軸、z軸軸方向に分けて3D表示します。109をかけて単位をnAとします。

rule4:=r[z]=0.01://1cm上
range:=(r[x]=0..0.2,r[y]=-0.1..0.1):
E2:=subs(E1[1]*1e9,rule4):
E3:=subs(E1[2]*1e9,rule4):
E4:=subs(E1[3]*1e9,rule4):
map([E2,E3,E4],plotfunc3d,range);

MuPAD graphics
MuPAD graphics
MuPAD graphics
[]

plot::VectorField3dを使ってベクトル場を計算させ、plot で表示させます。

range:=(r[x]=0.09..0.11,r[y]=-0.01..0.01,r[z]=-0.01..0.01)://系列
field:=plot::VectorField3d([E1[1],E1[2],E1[3]],range):
p1:=plot(field);

MuPAD graphics

 

電流の位置を(7,0,0)cm、大きさを(0,10,20)A、センサー位置を中心から半径12cmの球面上とし、球面上の法線方向の磁場の大きさを求めます。

rule1:=[Symbol::mu[0]=4*PI*1E-7]:
rule2:=[a[x]=r[x]-7E-2,a[y]=r[y],a[z]=r[z]]:
rule3:=[q[x]=0,q[y]=10,q[z]=20]:
R:=matrix([r[x],r[y],r[z]]):
Z:=R/norm(R,2):
Bz:=linalg::scalarProduct(B,Z,Real):
E1:=subs(Bz,[op(rule1),op(rule2),op(rule3)]):
E2:=simplify(E1);//まとめます。

-(0.00000001*(14.0*r[y] - 7.0*r[z]))/((abs(r[x])^2 + abs(r[y])^2 + abs(r[z])^2)^(1/2)*(abs(r[y])^2 + abs(r[z])^2 + abs(r[x] - 0.07)^2)^(3/2))

(rx,ry,rz)球座標に変換してplot::Surfaceを使います。
vertex上でBzの強度によって色を変える方法が分からないので球面上の凹凸で対応します。

rule4:=r[x]=0.12*cos(Symbol::theta)*cos(Symbol::phi):
rule5:=r[y]=0.12*sin(Symbol::theta)*cos(Symbol::phi):
rule6:=r[z]=0.12*sin(Symbol::phi):
E3:=(E2*R)*2000+R://補正
E4:=subs(E3,[rule4,rule5,rule6]):
range:=Symbol::theta=-PI..PI,Symbol::phi=-PI/2..PI/2:
p1:=plot::Surface(E4,range):
plot(p1);

MuPAD graphics