xbr.m

function BR=xbr(Rp,Zp,Rc,Zc,Sc)
%
% Compute radial field
%
% function BR=xbr(Rp,Zp,Rc,Zc,Sc);
% compute BR(T) at the points (Rp,Zp) due to unit currents (A) in (Rc,Zc) with area Sc
% BR(length(Rp),length(Rc)): each row refers to a point, each column to a coil
%       
RR   = Rp*ones(size(Rc')); 
ZZ   = Zp*ones(size(Zc')); 
rr   = ones(size(Rp))*Rc';
zz   = ones(size(Rp))*Zc';
ss   = ones(size(Rp))*Sc';
RA   = RR+rr;					% RAij= Ri+rj
DR   = RR-rr;					% DRij= Ri-rj
DZ   = ZZ-zz;					% DZij= Zi-zj
D2   = DR.*DR+DZ.*DZ+eps;
R2Z  = RR.*RR+DZ.*DZ;
K2   = 4*(Rp*Rc')./(RA.*RA+DZ.*DZ);
KI   = sqrt(RA.*RA+DZ.*DZ);
[K,E]= ellipke(K2);
BR   =-0.2e-6*(DZ./RR).*(K-(R2Z+rr.*rr)./D2.*E)./KI;

BR(RR==0) = 0;

inear = find(D2<ss/pi & D2 ~= eps);
if isempty(inear)==0;
  disp('some points near to coils!');
  BR(inear) = 2e-7*pi*DZ(inear)./ss(inear);
end

ising = find(D2 == eps);
if isempty(ising)==0;
  disp('+++ WARNING: some point overlaps with coil centre +++');
  BR(ising) = 0;
end
%
%
return