Question: How to calculate this hard integral?

How to calculate the integral of (z-z0)*z/sqrt((x-x0)^2+(y-y0)^2+(z-z0)^2)
over the unit sphere {(x,y,z):x^2+y^2+z^2<=1}
under the assumtion x0^2+y0^2+z0^2<=1 (x0^2+y0^2+z0^2>1)?
Its physical interpretation suggests the integral can be expressed through  elementary functions of the parameters.

My tries are
VectorCalculus:-int((z-z0)*z/sqrt((x-x0)^2+(y-y0)^2+(z-z0)^2),[x,y,z]=Sphere(<0,0,0>,1)) assuming x0^2+y0^2+z0^2<=1;

and

VectorCalculus:-int(eval((z-z0)*z/sqrt((x-x0)^2+(y-y0)^2+(z-z0)^2),
[x=r*sin(psi)*cos(theta),y=r*cos(psi)*sin(theta),z=r*cos(psi)])*r^2*sin(psi),
[r,psi,theta]=Parallelepiped(0..1,0..Pi,0..2*Pi)) assuming x0^2+y0^2+z0^2<=1;

The both are spinning on my comp. Also

VectorCalculus:-int((z-1/4)*z/sqrt((x-1/2)^2+(y-1/3)^2+(z-1/4)^2),[x,y,z]=Sphere(<0,0,0>,1),numeric);

is spinning.
Edt. The omitted part of the code assuming x0^2+y0^2+z0^2<=1 is added.