function B = infdipole(r,r0,muf)
B = [0 0 0];
zz = r-r0;
dd = (3.*zz*dot(zz,muf)/norm(zz)^5 - muf/norm(zz)^3);
B = 1/4/pi*dd;

return