stacksize(50000000);
loadmatfile('.../Dewar.mat');
getf('.../Sarvas.sce');
T=VectorView.triangle;P=VectorView.nodes;
//calculation flux by Sarvas's rule
q=[0.07,0,0,0,0,20*1e-9];Center=[0,0,40]/1000;
B=Sarvas(q,P',Center)*1e+15; // fT/cm
各三角メッシュ上での法線ベクトルは、頂点を有する三角平面の法線ベクトルとなりますが、 それだと磁場分布が凸凹になります。 そこで各三角平面の法線ベクトルを外積の大きさで重み付けすることにします。
// calculation of normal vector at each verteces
VP=zeros(P);
for x=1:size(P,2);...
k=(T(1,:)-x).*(T(2,:)-x).*(T(3,:)-x);...
k(k~=0)=1;k=(1:length(k)).*k;...
k(k==0)=[];TT=T;TT(:,k)=[];...
for y=1:size(TT,2);...
if TT(2,y)==x;TT(2,y)=TT(3,y);TT(3,y)=TT(1,y);TT(1,y)=x;...
elseif TT(3,y)==x;TT(3,y)=TT(2,y);TT(2,y)=TT(1,y);TT(1,y)=x;...
end;end;...
V1=P(:,TT(2,:))-P(:,TT(1,:));V2=P(:,TT(3,:))-P(:,TT(1,:));...
VP(:,x)=sum([V1(2,:).*V2(3,:)-V1(3,:).*V2(2,:);...
V1(3,:).*V2(1,:)-V1(1,:).*V2(3,:);...
V1(1,:).*V2(2,:)-V1(2,:).*V2(1,:)],2);end;
VP=VP./(ones(3,1)*sum(VP.^2,1).^0.5);VP=VP';
//drawing process
Bz=sum(B.*VP,2)';
ntri=size(T,2);
T=T($:-1:1,:); //裏返し
X=matrix(P(1,T),3,ntri);
Y=matrix(P(2,T),3,ntri);
Z=matrix(P(3,T),3,ntri);
Bz=matrix(Bz(T),3,ntri);
scale=32;scale2=scale*2;
Bz=Bz/250*scale+scale;Bz(Bz<1)=1;Bz(Bz>scale2)=scale2;
scf();set(gcf(),'color_map',jetcolormap(scale2));
set(gcf(),'visible','off');
plot3d(X,Y,Z);ha=gca();f=ha.children.data;
TL=tlist(['3d','x','y','z','color'],f.x,f.y,f.z,Bz);
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';ha.rotation_angles=[60,0];
set(gcf(),'visible','on');
まだギザギザがあります。 曲面fittingなどの技術があればもう少しsmoothに補間できるのかもしれません。
phong shadingの真似をします。
VP=VP';
Tn=zeros(size(VP,2),size(VP,2));
for x=1:size(T,2);...
TT=T(:,x);TT=TT(:)';...
Tn(TT(:),TT(:))=Tn(TT(:),TT(:))+1;...
end;
VPI=VP;
for x=1:size(VP,2);...
k=Tn(x,:);k(k>0)=1;k=(1:length(k)).*k;k(k==0)=[];...
VPI(:,x)=mean(VP(:,k),2);...
end;VPI=VPI';
Bz=sum(B.*VPI,2)';
Bz=matrix(Bz(T),3,ntri);
scale=32;scale2=scale*2;
Bz=Bz/250*scale+scale;Bz(Bz<1)=1;Bz(Bz>scale2)=scale2;
set(gcf(),'visible','off');
TL=tlist(['3d','x','y','z','color'],f.x,f.y,f.z,Bz);
ha.children.data=TL;
set(gcf(),'visible','on');
ほとんど変わりありません。残念。