function B=Bfield(x,y,z)
mu0=4*pi*10^-7;
rhat=[x y z]/(x^2+y^2+z^2)^.5;
m=[0 0 1];
%B=mu0/(x^2+y^2+z^2)^(3/2)/4/pi*(3*(m*rhat')*rhat-m);
B=[0,0,1];