vb_sarvas_sensor

 Magnetic field at MEG sensors for one dipole 
     in spherical brain model (Sarvas Eq.)
 B = vb_sarvas_sensor(P, Q, R, S)

 INPUT
   Dipole current source
   P : current dipole position   ( 1 x 3 )
   Q : current dipole moment     ( 1 x 3 )
   MEG sensor
   R : sensor position       ( NS x 3 )
   S : sensor orientation    ( NS x 3 )
 OUTPUT
   B : Lead field matrix         ( NS x 1 )

   NS : # of sensor
 
  2006-12-15 Made by M. Sato

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

0001 function    B = vb_sarvas_sensor(P, Q, R, S)
0002 % Magnetic field at MEG sensors for one dipole
0003 %     in spherical brain model (Sarvas Eq.)
0004 % B = vb_sarvas_sensor(P, Q, R, S)
0005 %
0006 % INPUT
0007 %   Dipole current source
0008 %   P : current dipole position   ( 1 x 3 )
0009 %   Q : current dipole moment     ( 1 x 3 )
0010 %   MEG sensor
0011 %   R : sensor position       ( NS x 3 )
0012 %   S : sensor orientation    ( NS x 3 )
0013 % OUTPUT
0014 %   B : Lead field matrix         ( NS x 1 )
0015 %
0016 %   NS : # of sensor
0017 %
0018 %  2006-12-15 Made by M. Sato
0019 %
0020 % Copyright (C) 2011, ATR All Rights Reserved.
0021 % License : New BSD License(see VBMEG_LICENSE.txt)
0022 
0023 % Dipole position vector
0024 P1 = P(1);
0025 P2 = P(2);
0026 P3 = P(3);
0027 
0028 % Dipole moment vector
0029 Q1 = Q(1);
0030 Q2 = Q(2);
0031 Q3 = Q(3);
0032 
0033 % Sensor position vector
0034 R1 = R(:,1);
0035 R2 = R(:,2);
0036 R3 = R(:,3);
0037 
0038 % Difference from dipole to sensor : (NS x 3)
0039 dP = [(R1-P1), (R2-P2), (R3-P3)];
0040 
0041 % Square norm    :  (NS x 1)
0042 dd = sum(dP.^2, 2);
0043 d  = sqrt(dd);
0044 rr = sum(R.^2,2);
0045 r  = sqrt(rr);
0046 
0047 % Inner product  : (NS x 3) * ( 3 x 1 ) = (NS x 1)
0048 dr = sum( dP .* R , 2);
0049 pr = R * P(:);
0050 
0051 % Cross product (Q x P) : (3 x 1)
0052 QP = [Q2.*P3-Q3.*P2; ...
0053       Q3.*P1-Q1.*P3; ...
0054       Q1.*P2-Q2.*P1  ];
0055 
0056 % Denominator of magnetic potencial  (NS x 1)
0057 f  = dd .* r + d .* dr ;
0058 
0059 % ( Gradient of 'f' ) * S
0060 df = (d + dd./r).* sum(R.*S,2)  + ( d + 2*r + dr./d ).* sum(dP.*S,2);
0061 
0062 % Magnetic field for sensor orientation vector 'S'
0063 B  = (S*QP)./f - (R*QP).*df./(f.^2);
0064 B  = (10^-7) * B;
0065 
0066 return
0067 % END of program
0068 
0069 %% Execution speed comparison with repmat
0070 %  Time
0071 %  0.23   dP = [(R1-P1), (R2-P2), (R3-P3)];
0072 %
0073 %% This code is faster than the following repmat code for large data
0074 %% since no memory allocation is needed
0075 % ( R1, R2, R3 are scalars)
0076 %
0077 %  0.36   R  = repmat(R',[size(P,1) 1]);
0078 %  0.17   dP = R - P;