% Forward calculation for MEG (homogenous sphere)

function ret = sphere_meg( r, u, 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.                                                                      

FACT = 1.0/3.18;
l = [2 3 1 2];  
    
ret(1)= 0.0;
ret(2)= 0.0;
ret(3)= 0.0;

rd(1)= r(1) - d(1);
rd(2)= r(2) - d(2);
rd(3)= r(3) - d(3);

rb= r(1)*r(1) + r(2)*r(2) + r(3)*r(3);
db= d(1)*d(1) + d(2)*d(2) + d(3)*d(3);
rdb= rd(1)*rd(1) + rd(2)*rd(2) + rd(3)*rd(3);
if( rb==0.0 | rdb==0.0 )
   return;
end;

rb= sqrt( rb );
rdb= sqrt( rdb );
if( db~= 0.0 )
	db= sqrt( db );
end;
   

alpha= -1.0 / rdb / ( rb*rdb + r(1)*rd(1)+r(2)*rd(2)+r(3)*rd(3) );
a= 1.0/rdb - 2.0*alpha*db*db - 1.0/rb;
b= 2.0 * alpha * ( r(1)*d(1)+r(2)*d(2)+r(3)*d(3) );
c= -1.0 / ( rdb*rdb*rdb ) * ( rd(1)*d(1)+rd(2)*d(2)+rd(3)*d(3) );
beta = a * ( r(1)*u(1)+r(2)*u(2)+r(3)*u(3) );
beta= beta + b * ( d(1)*u(1)+d(2)*u(2)+d(3)*u(3) );
beta= beta + c * ( rd(1)*u(1)+rd(2)*u(2)+rd(3)*u(3) );
rds= r(1)*d(1) + r(2)*d(2) + r(3)*d(3);
c= rb*rb * db*db - rds*rds;
if( c< 1.0e-10 )
	beta= 0.0;
else
   beta= beta/c;
end;

for i=1:3,   
   a= u( l(i) ) * d( l(i+1) ) - u( l(i+1) ) * d( l(i) );
   b= r( l(i) ) * d( l(i+1) ) - r( l(i+1) ) * d( l(i) );
   ret(i)= FACT * ( alpha * a + beta * b );
end;

return;

