傾斜磁場

Neuromag社の脳磁計は傾斜磁場を捉えています。 Sarvasの式を用い、baseline 28mm の傾斜磁場のsimulationをしてみます。
getf('Sarvas.sce');
stacksize(50000000);
t=linspace(-%pi,%pi,37)';
p=linspace(-%pi/2,%pi/2,19);
x=cos(t)*cos(p);
y=sin(t)*cos(p);
z=ones(t)*sin(p);
q=[0.07,0,0,0,0,20*1e-9];
r=0.12;L=[x(:),y(:),z(:)];
B=Sarvas(q,r*L,[0,0,0]);
B=sum(B.*L,2);
B=matrix(B,size(x,1),size(x,2));
[xx,yy,zz]=nf3d(x,y,z);
[BB,yy,zz]=nf3d(B,y,z);
xx=xx($:-1:1,:);yy=yy($:-1:1,:);zz=zz($:-1:1,:);BB=BB($:-1:1,:);
scale=64/2;scale2=scale*2;
BB=(BB*1e+15)/250*scale+scale;
BB(BB<1)=1;BB(BB>scale2)=scale2;
scf();set(gcf(),'color_map',jetcolormap(scale*2));
set(gcf(),'visible','off');
plot3d(xx,yy,zz);
ha=gca();f=ha.children.data;
TL=tlist(['3d','x','y','z','color'],f.x,f.y,f.z,BB);
ha.children.data=TL;
ha.children.color_mode=-1;
ha.children.color_flag=3;
ha.isoview='on';ha.box='off';ha.margins=[0,0,0,0];
ha.axes_visible='off';set(gcf(),'visible','on');

球面上の点を

とおくとき、(x,y,z)における接平面の 横方向θの向きは、

同じく接平面の縦方向φ方向の向きは

となります。
横方向の傾斜磁場を求めます。
dx=-sin(t)*cos(p);
dy=cos(t)*cos(p);
dz=zeros(dx);
dR=(dx.^2+dy.^2+dz.^2).^0.5;
dx=dx./dR;dy=dy./dR;dz=dz./dR;
dR=[dx(:),dy(:),dz(:)];
wd=14/1000;
Bs1=Sarvas(q,r*L+wd*dR,[0,0,0]);
Bs2=Sarvas(q,r*L-wd*dR,[0,0,0]);
dBx=matrix(sum((Bs2-Bs1).*L,2),size(x,1),size(x,2));
dBx=dBx/2/wd;
dBxc=dBx/(1e-11)*scale+scale;
dBxc(dBxc<1)=1;dBxc(dBxc>scale2)=scale2;
[BB,yy,zz]=nf3d(dBxc,y,z);
yy=yy($:-1:1,:);zz=zz($:-1:1,:);BB=BB($:-1:1,:);
scf();set(gcf(),'color_map',jetcolormap(scale*2));
set(gcf(),'visible','off');
plot3d(xx,yy,zz);
ha=gca();f=ha.children.data;
TL=tlist(['3d','x','y','z','color'],f.x,f.y,f.z,BB);
ha.children.data=TL;
ha.children.color_mode=-1;
ha.children.color_flag=3;
ha.isoview='on';ha.box='off';ha.margins=[0,0,0,0];
ha.axes_visible='off';set(gcf(),'visible','on');

縦方向の傾斜磁場を求めます。
dx=-cos(t)*sin(p);
dy=-sin(t)*sin(p);
dz=ones(t)*cos(p);
dR=(dx.^2+dy.^2+dz.^2).^0.5;
dx=dx./dR;dy=dy./dR;dz=dz./dR;
dR=[dx(:),dy(:),dz(:)];
wd=14/1000;
Bs1=Sarvas(q,r*L+wd*dR,[0,0,0]);
Bs2=Sarvas(q,r*L-wd*dR,[0,0,0]);
dBy=matrix(sum((Bs2-Bs1).*L,2),size(x,1),size(x,2));
dBy=dBy/2/wd;
dBxc=dBy/(1e-11)*scale+scale;
dBxc(dBxc<1)=1;dBxc(dBxc>scale2)=scale2;
[BB,yy,zz]=nf3d(dBxc,y,z);
yy=yy($:-1:1,:);zz=zz($:-1:1,:);BB=BB($:-1:1,:);
scf();set(gcf(),'color_map',jetcolormap(scale*2));
set(gcf(),'visible','off');
plot3d(xx,yy,zz);
ha=gca();f=ha.children.data;
TL=tlist(['3d','x','y','z','color'],f.x,f.y,f.z,BB);
ha.children.data=TL;
ha.children.color_mode=-1;
ha.children.color_flag=3;
ha.isoview='on';ha.box='off';ha.margins=[0,0,0,0];
ha.axes_visible='off';set(gcf(),'visible','on');

θ方向とφ方向の傾斜磁場の二乗和根を求めます。
BB=(dBx.^2+dBy.^2).^0.5;
BB=BB/(1e-11)*scale+scale;
BB(BB<1)=1;BB(BB>scale2)=scale2;
[BB,yy,zz]=nf3d(BB,y,z);
yy=yy($:-1:1,:);zz=zz($:-1:1,:);BB=BB($:-1:1,:);
scf();set(gcf(),'color_map',jetcolormap(scale*2));
set(gcf(),'visible','off');
plot3d(xx,yy,zz);
ha=gca();f=ha.children.data;
TL=tlist(['3d','x','y','z','color'],f.x,f.y,f.z,BB);
ha.children.data=TL;
ha.children.color_mode=-1;
ha.children.color_flag=3;
ha.isoview='on';ha.box='off';ha.margins=[0,0,0,0];
ha.axes_visible='off';set(gcf(),'visible','on');

NeuromagのSource ModellingのGradient Isofield Contour mapは傾斜磁場の二乗根和です。 マイナスがないので、等磁場線はプラスの赤だけとなります。
θ方向とφ方向で２回偏微分した場合は

となります。縦方向・横方向で偏微分した場合の傾斜磁場を求めます。
dx=sin(t)*sin(p);
dy=-cos(t)*sin(p);
dz=0;
dR=(dx.^2+dy.^2+dz.^2).^0.5;
dR(dR==0)=1;
dx=dx./dR;dy=dy./dR;dz=dz./dR;
dR=[dx(:),dy(:),dz(:)];
wd=14/1000;
Bs1=Sarvas(q,r*L+wd*dR,[0,0,0]);
Bs2=Sarvas(q,r*L-wd*dR,[0,0,0]);
dBy=matrix(sum((Bs2-Bs1).*L,2),size(x,1),size(x,2));
dBy=dBy/2/wd;
dBxc=dBy/(1e-11)*scale+scale;
dBxc(dBxc<1)=1;dBxc(dBxc>scale2)=scale2;
[BB,yy,zz]=nf3d(dBxc,y,z);
yy=yy($:-1:1,:);zz=zz($:-1:1,:);BB=BB($:-1:1,:);
scf();set(gcf(),'color_map',jetcolormap(scale*2));
set(gcf(),'visible','off');
plot3d(xx,yy,zz);
ha=gca();f=ha.children.data;
TL=tlist(['3d','x','y','z','color'],f.x,f.y,f.z,BB);
ha.children.data=TL;
ha.children.color_mode=-1;
ha.children.color_flag=3;
ha.isoview='on';ha.box='off';ha.margins=[0,0,0,0];
ha.axes_visible='off';set(gcf(),'visible','on');