function Bfinal=ChampSphereLT(x,y,theta,phi,a,d,M)

mu_0 = 4*pi*10^-7;

%%->Magnetic field on a plane from a magnetized sphere 
%%->The direction of magnetization is imposed by the external magnetic field B0
%%->(x,y,z) is defined by the diamond crystal structure axes
%%->a=Sphere Radius (in m)
%%->d=distance of sphere edge from NV plane
%%->M Magnetization

P=[x;y;a+d]; %Point (x,y) on the NV centers plane

%We want to express the coordinates of P in the frame defined by the
%Magnetization axis -> we need a rotation from the frame of the crystal to
%the frame with z'//M -> Rotation of angle phi (angle between B0 and z)
%around the unit vector perpendicular to B0xy

%%Definition of the rotation matrix
u=[-cos(theta);sin(theta);0];%Definition of the rotation axis u=VectorProduct(B0,z) 
c=cos(-phi);
s=sin(-phi);
ux=u(1);
uy=u(2);
uz=u(3);
R=[ux^2*(1-c)+c,ux*uy*(1-c)-uz*s,ux*uz*(1-c)+uy*s;...
    ux*uy*(1-c)+uz*s,uy^2*(1-c)+c,uy*uz*(1-c)-ux*s;...
    ux*uz*(1-c)-uy*s,uy*uz*(1-c)+ux*s,uz^2*(1-c)+c]; %Cf Wikipedia for rotation of angle theta along u
%%

Pprime=R*P;%Point Coordinate in the rotated frame

%%Conversion cartesian coordinates -> spherical coordinates
rPprime=norm(Pprime);%r
thetaPprime=acos(Pprime(3)/norm(Pprime));%theta
if thetaPprime==0,phiPprime=0;%Singular case
elseif Pprime(2)>=0
phiPprime=acos(Pprime(1)/(rPprime*sin(thetaPprime)));%Phi
else 
    phiPprime=-acos(Pprime(1)/(rPprime*sin(thetaPprime)));%Phi
end

% %%Conversion cartesian coordinates -> spherical coordinates
% rPprime=norm(Pprime);%r
% phiPprime=acos(Pprime(3)/norm(Pprime));%theta
% if phiPprime==0,thetaPprime=0;%Singular case
% elseif Pprime(2)>=0
%     thetaPprime=acos(Pprime(1)/(rPprime*sin(phiPprime)));%Phi
% else 
%     thetaPprime=-acos(Pprime(1)/(rPprime*sin(phiPprime)));%Phi
% end

%%Calculation of B', in the rotated frame (cf Jackson : in this the magnetization is parallel to z' by definition)
BrPrime=mu_0*M*2/3*(a/rPprime).^3*cos(thetaPprime);
BthetaPrime=mu_0*M*1/3*(a/rPprime).^3*sin(thetaPprime);

%%Back Conversion spherical -> cartesian
BxPrime=(BrPrime*sin(thetaPprime)+BthetaPrime*cos(thetaPprime))*cos(phiPprime);
ByPrime=(BrPrime*sin(thetaPprime)+BthetaPrime*cos(thetaPprime))*sin(phiPprime);
BzPrime=cos(thetaPprime)*BrPrime-sin(thetaPprime)*BthetaPrime;

%%Inverse rotation to go back to the diamond frame
Bfinal=R\[BxPrime;ByPrime;BzPrime]*1000;
