function dDelta_B=Delta_B(dXo,dYo,dZo,... % Observer coordinate
dA,... % X azimuth angle coordinate
dXc,dYc,dZc,... % Sphere coordinate
dR,... % Sphere radius
dK,... % Magnetic susceptibility
dMr,dIr,dDr,... % Remanence sphere parameters
dB0,dI,dD) % Background field parameters
MC=4*pi*1E-7;% permeability of vacuum
DEG2ARC=pi/180.0;
NT2T=1E-9;
dA=dA*DEG2ARC;
dIr=dIr*DEG2ARC;
dDr=dDr*DEG2ARC;
dI=dI*DEG2ARC;
dD=dD*DEG2ARC;
dB0=dB0*NT2T;
dLx=cos(dIr)*cos(dDr-dA);
dLy=cos(dIr)*sin(dDr-dA);
dLz=sin(dIr);
dNx=cos(dI)*cos(dD-dA);
dNy=cos(dI)*sin(dD-dA);
dNz=sin(dI);
% Total magnetization three component
dMi=dK*dB0/MC; %Size of induced magnetization
dMx=dMi*dNx+dMr*dLx;
dMy=dMi*dNy+dMr*dLy;
dMz=dMi*dNz+dMr*dLz;
%Two partial derivatives of the gravitational potential of Mathematics
dVxx=Vaa(dXo,dYo,dZo,dXc,dYc,dZc);
dVyy=Vaa(dYo,dZo,dXo,dYc,dZc,dXc);
dVzz=Vaa(dZo,dXo,dYo,dZc,dXc,dYc);
dVxy=Vab(dXo,dYo,dZo,dXc,dYc,dZc);
dVyz=Vab(dYo,dZo,dXo,dYc,dZc,dXc);
dVzx=Vab(dZo,dXo,dYo,dZc,dXc,dYc);
dV=4*pi*dR^3/3.0;
%Three component of magnetic anomaly
dBax=MC*dV*(dMx*dVxx+dMy*dVxy+dMz*dVzx)/(4*pi);
dBay=MC*dV*(dMx*dVxy+dMy*dVyy+dMz*dVyz)/(4*pi);
dBaz=MC*dV*(dMx*dVzx+dMy*dVyz+dMz*dVzz)/(4*pi);
dDelta_B=dBax*dNx+dBay*dNy+dBaz*dNz;