function lead=lf1(e,d);
%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      
% computes the leadfield relating sensor e with solution point d
% assumes conductivity=1, there is no factor 1/4*pi*sigma, sphere with radio=1.
% e=(xe,ye,ze) electrode coordinates.
% d=(xd,yd,zd) dipole or solution points coordinates.
% lead=(lx,ly,lz) resulting leadfield components.

% reprogramed 1-1-97 from Delia version to be faster 
% is really faster in Matlab?  not checked!

e=e(:); d=d(:);  %column vectors

dn=norm(e-d,2);
d3=2/dn^3;
ol=sum(e.*d);
dd1s=dn*(dn + 1.0 - ol);

factore=d3+ (1/dn) + (ol/dd1s);
factord=d3+ (1/dd1s);

lead=(factore*e-factord*d);

%over
