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.