function R = calBdB(P)
%
% INPUT  : P = coordinates of the voxel in the slice frame
% OUTPUT : 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) = proton density in the domain P belongs to.
%
% (c) Stephane  Balac - MAPLY - INSA Lyon - mars 2001
%
global r1    % parameters of the experiment
global r2
global mu
global mu0
global rho
global CA
global CB
global CC
global CD
global nr
global np
global adim

xr = P(1);
xp = P(2);
r = sqrt(xr^2 + xp^2);
if  (r < r1)
% P belongs to domain 1
  R(1) = - mu0*mu(1)*CA;
  R(2) = calr2(-mu0*mu(1)*CA,0,nr(1),np(1),xp,xr);
  R(3) = rho(1);
elseif (r < r2) 
% P belongs to domain 2
  R(1) = - mu0*mu(2)*(CC - CD*(2*P(1)^2/r^4 - 1/r^2));
  R(2) =  calr2(-mu0*mu(2)*CC,mu0*mu(2)*CD,nr(1),np(1),xp,xr);
  R(3) = rho(2);
else
  R(1) =  mu0*mu(3)*CB*(2*P(1)^2/r^4 - 1/r^2);
  R(2) = calr2(0,mu0*mu(3)*CB,nr(1),np(1),xp,xr);
  R(3) = rho(3);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function R2 = calr2(alpha1,alpha2,ax,bx,xp,xr)
numer=-2*alpha2*(2*ax^2-1)*xr^3-12*alpha2*ax*bx*xp*xr^2 ...
-2*alpha2*(4*bx^2-2*ax^2-1)*xp^2*xr+4*alpha2*ax*bx*xp^3;
denom=(xr^2+xp^2)^3;
R2 = numer/denom;

