function B=Sarvas(Q,S,p);
//Q = (Nq,[x,y,z,qx,qy,qz]);//dipole position & moment
//S = (Ns,[x,y,z]);         //Sensor position
//p = [x,y,z];              //sphere center
Ns=size(S,1);
if size(Q,1)==1;
    Q=ones(Ns,1)*Q;
end;
P=ones(Ns,1)*p;
R=S-P;
R0=Q(:,1:3)-P;
A=R-R0;
r=sum(R.^2,2).^0.5;
r0=sum(R0.^2,2).^0.5;
a=sum(A.^2,2).^0.5;
ar=sum(A.*R,2);
F=a.*(r.*a+ar);
n=a+2*r+ar./a;
nn=(a.^2)./r+a+n;
n=n*ones(1,3);
nn=nn*ones(1,3);
NablaF=nn.*R-n.*R0;//vector!
Q=Q(:,4:6);
QR0=zeros(Q);
QR0(:,1)=Q(:,2).*R0(:,3)-Q(:,3).*R0(:,2);
QR0(:,2)=Q(:,3).*R0(:,1)-Q(:,1).*R0(:,3);
QR0(:,3)=Q(:,1).*R0(:,2)-Q(:,2).*R0(:,1);
B=zeros(Q);
F=F*ones(1,3);
QR0R=sum(QR0.*R,2)*ones(1,3);//vector!
B=QR0-(QR0R.*NablaF)./F;
B=B./F*1e-7;
endfunction;


