Sarvasの式2

Sarvasの式で、球の中心を動かしてみます。
電流の位置を(7,0,0)cm、大きさを(0,0,20)nAとして中心(0,0,0)cm、半径12cmの球面上のセンサーでの磁場を計算しますが、 均質導体モデルの球の中心座標をいろいろ変化させてみます。

getf('../Sarvas.sce');
stacksize(50000000);
t=linspace(-%pi,%pi,37)';
p=linspace(-%pi/2,%pi/2,19);
r=0.12;
x=r*cos(t)*cos(p);
y=r*sin(t)*cos(p);
z=r*ones(t)*sin(p);
scf();
q=[0.07,0,0,0,0,20*1e-9];
k=prod(size(x));
x=matrix(x,k,1);
y=matrix(y,k,1);
z=matrix(z,k,1);
R=[x,y,z];

//// ルーチンタスクの関数化 ////
function ShowFlux(R,q,B);
set(gcf(),'visible','off');
param3d(R(:,1),R(:,2),R(:,3));
h=gca();h.box='off';h.margins=[0,0,0,0];
h.children.mark_mode='on';
h.children.mark_size=3;
Q=ones(2,1)*q(1:3);
param3d(Q(:,1),Q(:,2),Q(:,3));
h.children(1).foreground=5;
h.children(1).mark_mode='on';
h.children(1).mark_size=3;
Q(2,:)=q(1:3)+q(4:6)*4e+6;
param3d(Q(:,1),Q(:,2),Q(:,3));
h.children(1).foreground=5;
B=B*1e+11;
Fx=[R(:,1),R(:,1)+B(:,1)]';
Fy=[R(:,2),R(:,2)+B(:,2)]';
Fz=[R(:,3),R(:,3)+B(:,3)]';
set(gcf(),'visible','off');
k=prod(size(R(:,1)));
for x=1:k;param3d(Fx(:,x),Fy(:,x),Fz(:,x));end;
h.isoview='on';
set(gcf(),'visible','on');
endfunction;

//// Draw Flux Process //// 結構時間がかかります
subplot(221); ShowFlux(R,q,Sarvas(q,R,[0,0,0])); // (0,0,0)cm
subplot(222); ShowFlux(R,q,Sarvas(q,R,[0.04,0,0])); // (4,0,0)cm
subplot(223); ShowFlux(R,q,Sarvas(q,R,[0,0.04,0])); // (0,4,0)cm
subplot(224); ShowFlux(R,q,Sarvas(q,R,[0,0,0.04])); // (0,0,4)cm


左上は(0,0,0)cmです。Biot-Sarvartと同じです。
右上は(4,0,0)cm、左下は(0,4,0)cm、右下は(0,0,4)cmです。
Biot-Savartの式では全く関係ありませんが、導体球モデルの位置により磁力線が変化することがわかりました。

新たなsubplotが追加されるたびに描画速度が落ちているような気がします。 subplot数を6個にすると絶望的な遅さになるので4個にとどめました。