Home > functions > leadfield > bem > vb_sensor_mag.m

vb_sensor_mag

PURPOSE ^

境界三角形 (x1, x2, x3) の単位電流 が X に作る磁場の Q 方向射影

SYNOPSIS ^

function BS = vb_sensor_mag(Vs,Fs,s0,X,Q)

DESCRIPTION ^

 境界三角形 (x1, x2, x3) の単位電流 が X に作る磁場の Q 方向射影
 Vs : 境界三角面の座標
 Fs : 境界三角面頂点インデックス
 s0 : 三角面の外向き面積法線
 X  : 磁気センサー座標             (Npick,  3)
 Q  : 磁場の観測方向               (Npick,  3)

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function    BS = vb_sensor_mag(Vs,Fs,s0,X,Q)
0002 % 境界三角形 (x1, x2, x3) の単位電流 が X に作る磁場の Q 方向射影
0003 % Vs : 境界三角面の座標
0004 % Fs : 境界三角面頂点インデックス
0005 % s0 : 三角面の外向き面積法線
0006 % X  : 磁気センサー座標             (Npick,  3)
0007 % Q  : 磁場の観測方向               (Npick,  3)
0008 %
0009 % Copyright (C) 2011, ATR All Rights Reserved.
0010 % License : New BSD License(see VBMEG_LICENSE.txt)
0011 
0012 % 2004-02-06 Taku Yoshioka
0013 % 2004-12-26 M. Sato modified
0014 
0015 %  変数の次元
0016 NS = size(X,1);        % 磁気センサー
0017 NV = size(Vs,1);    % 境界頂点
0018 NF = size(Fs,1);    % 境界三角面
0019 
0020 BS = zeros(NV,NS);  % 境界頂点上の単位電流が作る磁場
0021 
0022 % 三角面の頂点　triangle vertex
0023 x1 = Vs(Fs(:,1),:);
0024 x2 = Vs(Fs(:,2),:);
0025 x3 = Vs(Fs(:,3),:);
0026 
0027 %%% 境界頂点単位電流 が 観測点 に作る磁場計算 %%%
0028 
0029 for i=1:NS,
0030   % i-番目の観測点
0031   G0 = Q(i,:)';
0032   XP = X(i,:);
0033   XX = XP(ones(NF,1),:);
0034     
0035   % 境界三角形 (x1, x2, x3) 単位電流 が XP に作る磁場の G0 方向射影
0036   
0037   % センサ位置に対する三角面頂点の相対位置
0038   y1 = x1 - XX;
0039   y2 = x2 - XX;
0040   y3 = x3 - XX;
0041   
0042   % y1*((y2-y1)x(y3-y1)) が三角形の外向き法線になるように
0043   % y2 と y3 を入れ替え
0044   ix = find( sum(vb_cross2(y2-y1,y3-y1).*y1, 2) < 0);
0045   tmp = y2(ix,:);
0046   y2(ix,:) = y3(ix,:);
0047   y3(ix,:) = tmp;
0048   Fn = Fs; 
0049   Ftemp    = Fn(ix,2);
0050   Fn(ix,2) = Fn(ix,3);
0051   Fn(ix,3) = Ftemp;
0052 
0053   % 頂点相対位置ベクトルのノルム
0054   yy1 = sqrt(sum(y1.^2 ,2));
0055   yy2 = sqrt(sum(y2.^2 ,2));
0056   yy3 = sqrt(sum(y3.^2 ,2));
0057   
0058   % 双対基底ベクトル
0059   z1 = vb_cross2(y2,y3);
0060   z2 = vb_cross2(y3,y1);
0061   z3 = vb_cross2(y1,y2);
0062   
0063   d  = sum(z1.*y1,2);
0064   s  = z1+z2+z3;
0065   ss = sqrt(sum(s.^2,2));
0066   n  = s./ss(:,ones(1,3));
0067 
0068   % 三角形の辺
0069   w1 = y2 - y3;
0070   w2 = y3 - y1;
0071   w3 = y1 - y2;
0072 
0073   % 辺のノルム
0074   ww1 = sqrt(sum(w1.^2 ,2));
0075   ww2 = sqrt(sum(w2.^2 ,2));
0076   ww3 = sqrt(sum(w3.^2 ,2));
0077 
0078   wy1 = ( ww3.*yy2 - sum(w3.*y2, 2) )./( ww3.*yy1 - sum(w3.*y1, 2) );
0079   wy2 = ( ww1.*yy3 - sum(w1.*y3, 2) )./( ww1.*yy2 - sum(w1.*y2, 2) );
0080   wy3 = ( ww2.*yy1 - sum(w2.*y1, 2) )./( ww2.*yy3 - sum(w2.*y3, 2) );
0081 
0082   g1 = log(wy1)./ww3;
0083   g2 = log(wy2)./ww1;
0084   g3 = log(wy3)./ww2;
0085     
0086   cc = sum(y1.*n,2);
0087   c = cc(:,ones(1,3)).*n;
0088   c1 = y1-c;
0089   c2 = y2-c;
0090   c3 = y3-c;
0091 
0092   % 三角面の立体角
0093   omega = 2*atan(d./((yy1.*yy2.*yy3) ...
0094                     + yy1.*sum(y2.*y3,2) ...
0095                     + yy2.*sum(y3.*y1,2) ...
0096                     + yy3.*sum(y1.*y2,2)));
0097   
0098   g = sign(sum(s0.*n,2)) ...
0099               .*(g1.*sum(n.*vb_cross2(c1,c2),2) ...
0100                + g2.*sum(n.*vb_cross2(c2,c3),2) ...
0101                + g3.*sum(n.*vb_cross2(c3,c1),2) ...
0102                - sum(n.*c,2).*omega)./ss;
0103 
0104   % 各頂点の寄与
0105   omegaF1 = w1.*g(:,ones(1,3));
0106   omegaF2 = w2.*g(:,ones(1,3));
0107   omegaF3 = w3.*g(:,ones(1,3));
0108   
0109   for j = 1:NF
0110     BS(Fn(j,1),i) = BS(Fn(j,1),i)+omegaF1(j,:)*G0;
0111     BS(Fn(j,2),i) = BS(Fn(j,2),i)+omegaF2(j,:)*G0;
0112     BS(Fn(j,3),i) = BS(Fn(j,3),i)+omegaF3(j,:)*G0;
0113   end
0114 end;
0115 
0116 % BS(i,k) : 境界点 V(k,:) の単位電流 が X(i,:) に作る磁場
0117 %         : (Nsensor x Nvertex)
0118 BS = BS';
Generated on Tue 27-Aug-2013 11:46:04 by m2html © 2005