function R = calBdB(P)
% IN  : P = coordinates of the voxel
% OUT : the following values
%  R(1) = component along Z of the magnetic field
%             induced by the ball B'_z in P
%  R(2) = derivative along n_r of B'_z in P
%             d[B'_z(P)]/d[x_r]
%  R(3) = sub-domain number
%           0 if P belong to the ball
%           1 otherwise
global r0
global beta
global c
global Z

r = sqrt( P(1)^2 + P(2)^2 + (P(3)+c)^2 );
zr   = dot(Z,P) + Z(3)*c;

R(1) = (beta/r^5) * (3*zr^2 - r^2);
R(2) = (beta/r^7) * (6*Z(1)*zr*r^2 - 15*P(1)*zr^2 + 3*P(1)*r^2);
R(3) = 1; if r < r0, R(3) = 0; end      
