vb_dipole_field_magdipole

 Lead field matrix for magnetic dipole 
       (Note: it is different from current dipole)
   B = vb_dipole_field_magdipole(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

 Magnetic field produced by magnetic dipole
 Coil current can be approximated by magnetic dipole.

 磁気双極子磁場 (コイルに流れる電流の近似式)
 x0 の磁気双極子 J0 が 観測点 X に作る磁場の Q 方向射影

 B = 3(Q*Xd)(J0*Xd)/|Xd|^5 - (Q*J0)/|Xd|^3 , Xd = X - x0

  Ver-1.0 2007-7-9 Made by M. Sato

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

0001 function    B = vb_dipole_field_magdipole(x0,J0,X,Q)
0002 % Lead field matrix for magnetic dipole
0003 %       (Note: it is different from current dipole)
0004 %   B = vb_dipole_field_magdipole(x0,J0,X,Q)
0005 %
0006 % INPUT
0007 %   -Dipole current source
0008 %   x0 : current dipole position   ( NP x 3 )
0009 %   J0 : current dipole moment     ( NP x 3 )
0010 %   -MEG sensor
0011 %   X : sensor position            ( NS x 3 )
0012 %   Q : sensor orientation         ( NS x 3 )
0013 % OUTPUT
0014 %   B : Lead field matrix          ( NP x NS )
0015 %
0016 %   NP : # of dipole points
0017 %   NS : # of MEG sensor
0018 %
0019 % Magnetic field produced by magnetic dipole
0020 % Coil current can be approximated by magnetic dipole.
0021 %
0022 % 磁気双極子磁場 (コイルに流れる電流の近似式)
0023 % x0 の磁気双極子 J0 が 観測点 X に作る磁場の Q 方向射影
0024 %
0025 % B = 3(Q*Xd)(J0*Xd)/|Xd|^5 - (Q*J0)/|Xd|^3 , Xd = X - x0
0026 %
0027 %  Ver-1.0 2007-7-9 Made by M. Sato
0028 %
0029 % Copyright (C) 2011, ATR All Rights Reserved.
0030 % License : New BSD License(see VBMEG_LICENSE.txt)
0031 
0032 NP = size(x0,1);
0033 NS = size(X,1);
0034 
0035 if NP >= NS
0036     B  = zeros(NP,NS);
0037     for n=1:NS
0038         Xn = X(n,:);
0039         Qn = Q(n,:);
0040         Xd = vb_repadd( -x0, Xn);    % [NP x 1]
0041         RR = sum( Xd.^2, 2 );
0042 
0043         B(:,n) = (3*(Xd*Qn').*sum(J0.*Xd,2)./RR - (J0*Qn'))./(RR.^(3/2));
0044     end;
0045 else
0046     B  = zeros(NS,NP);
0047     for m=1:NP
0048         Jm = J0(m,:);
0049         xm = x0(m,:);
0050         Xd = vb_repadd( X, -xm);    % [NS x 1]
0051         RR = sum( Xd.^2, 2 );
0052 
0053         B(:,m) = (3*sum(Xd.*Q,2).*(Xd*Jm')./RR - (Q*Jm'))./(RR.^(3/2));
0054     end;
0055     B = B';
0056 end;
0057 
0058 B  = (10^-7) * B;