vb_dipole_magnetic

 Lead field matrix for dipole current (Biot-Savart)
   B = vb_dipole_magnetic(x0,J0,X,Q)

 INPUT
   Dipole current source
   x0 : current dipole position   ( NP x 3 )
   J0 : current dipole moment     ( NP x 3 )
   MEG sensor
   X : sensor position            ( NS x 3 )
   Q : sensor orientation         ( NS x 3 )
 OUTPUT
   B : Lead field matrix          ( NP x NS )

   NP : # of dipole points
   NS : # of MEG sensor

 Biot-Savart's law
 双極子磁場計算
 x0 の双極子 J0 が 観測点 X に作る磁場の Q 方向射影
 B = Q*(J0 x Xd)/|Xd|^3 , Xd = X - x0

  Ver-1.0 2004-12-22 Made by M. Sato

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

0001 function    B = vb_dipole_magnetic(x0,J0,X,Q)
0002 % Lead field matrix for dipole current (Biot-Savart)
0003 %   B = vb_dipole_magnetic(x0,J0,X,Q)
0004 %
0005 % INPUT
0006 %   Dipole current source
0007 %   x0 : current dipole position   ( NP x 3 )
0008 %   J0 : current dipole moment     ( NP x 3 )
0009 %   MEG sensor
0010 %   X : sensor position            ( NS x 3 )
0011 %   Q : sensor orientation         ( NS x 3 )
0012 % OUTPUT
0013 %   B : Lead field matrix          ( NP x NS )
0014 %
0015 %   NP : # of dipole points
0016 %   NS : # of MEG sensor
0017 %
0018 % Biot-Savart's law
0019 % 双極子磁場計算
0020 % x0 の双極子 J0 が 観測点 X に作る磁場の Q 方向射影
0021 % B = Q*(J0 x Xd)/|Xd|^3 , Xd = X - x0
0022 %
0023 %  Ver-1.0 2004-12-22 Made by M. Sato
0024 %
0025 % Copyright (C) 2011, ATR All Rights Reserved.
0026 % License : New BSD License(see VBMEG_LICENSE.txt)
0027 
0028 NP = size(x0,1);
0029 NS = size(X,1);
0030 
0031 B  = zeros(NP,NS);
0032 
0033 if NP >= NS
0034     for n=1:NS
0035         Xn = X(n,:);
0036         Qn = Q(n,:);
0037         QQ = Qn(ones(NP,1),:);
0038         Xd = Xn(ones(NP,1),:) - x0;
0039         Rx = sqrt(sum( Xd.^2, 2 )).^3;
0040             
0041         % 外積計算 JJ x Xd
0042         JB = [ J0(:,2).*Xd(:,3) - J0(:,3).*Xd(:,2), ...
0043                J0(:,3).*Xd(:,1) - J0(:,1).*Xd(:,3), ...
0044                J0(:,1).*Xd(:,2) - J0(:,2).*Xd(:,1)];
0045             
0046         B(:,n) = sum( JB.*QQ, 2 )./Rx;
0047     end;
0048 else
0049     for n=1:NP
0050         Jn = J0(n,:);
0051         xn = x0(n,:);
0052         JJ = Jn(ones(NS,1),:);
0053         Xd = X - xn(ones(NS,1),:);
0054         Rx = sqrt(sum( Xd.^2, 2 )).^3;
0055             
0056         % 外積計算 JJ x Xd
0057         JB = [ JJ(:,2).*Xd(:,3) - JJ(:,3).*Xd(:,2), ...
0058                JJ(:,3).*Xd(:,1) - JJ(:,1).*Xd(:,3), ...
0059                JJ(:,1).*Xd(:,2) - JJ(:,2).*Xd(:,1)];
0060             
0061         B(n,:) = (sum( JB.*Q, 2 )./Rx)';
0062     end;
0063 end;
0064 
0065 B  = (10^-7) * B;