Sarvasの式

ニューロマグ社の脳磁計では、Biot-Savartの式でなく、 均質導体球モデルであるSarvasの式が用いられています。

ここで B(r)は磁場を、 μ₀は透磁率(4π×10^-7N/A²)を、 Qは電流を、 r₀は球の中心に対する電流の位置ベクトルを、 rは球の中心に対するセンサーの位置ベクトルを、 aは電流の位置に対する、センサーの位置ベクトルを意味します。
Biot-Savartの式と比べ少し複雑になり、新たに球の中心座標の情報が必要になっています。
そこで以下のようなSCEファイル、Sarvas.sceを作成しました。

演算数を少なくするため、少し変形して

として計算しています。

になったらどうなるのだろうか？という疑問はありますが、無視しました。
シミュレーションしてみます。電流の位置を(7,0,0)cm、 大きさを(0,0,20)nA、球の中心を(0,0,0)cm、磁場計測地点を半径12cmの円周上とします。
t=linspace(-%pi,%pi,361)';
r=0.12;
x=r*cos(t);y=r*sin(t);z=zeros(x);
q=[0.07,0,0,0,0,20*1e-9];//[m] & [nA]
getf('BiotSavart.sce');
getf('Sarvas.sce');
Bb=BiotSavart(q,[x,y,z])*1e+15;//[fT]
Bs=Sarvas(q,[x,y,z],[0,0,0])*1e+15;//[fT]
scf(); subplot(121);plot2d(t/%pi*180,Bb);
h=gca();h.data_bounds=[-180,-400;180,800];
h.tight_limits='on';h.grid=[0,0];
h.tight_limits='on';h.grid=[0,0];
xtitle('Biot-Savart''s rule');
h.title.font_size=3;
subplot(122);plot2d(t/%pi*180,Bs);
h=gca();h.data_bounds=[-180,-400;180,800];
h.tight_limits='on';h.grid=[0,0];
xtitle('Sarvas''s rule')
h.title.font_size=3;

続いて法線方向の成分を計算します。
Bbz=sum(Bb.*[x,y,z],2);
Bsz=sum(Bs.*[x,y,z],2);
scf();
subplot(121);plot2d(t/%pi*180,Bbz);
h=gca();h.grid=[0,0];h.data_bounds=[-180,-30;180,30];
h.tight_limits='on';xtitle('Biot-Savart''s rule');
h.title.font_size=3;
subplot(122);plot2d(t/%pi*180,Bsz);
h=gca();h.grid=[0,0];h.data_bounds=[-180,-30;180,30];
h.tight_limits='on';xtitle('Sarvas''s rule');
h.title.font_size=3;

ん？
max(abs(Bbz-Bsz))=1.776D-14
まったく同じ結果になってしまいました。