I am trying to integrate the product of a function and a partially constant delta function over a sphere of constant radius $r$.
The integral is of the form
$$\int^{2\pi}_0 \int^{\pi}_0 f(\mathbf{r})\delta(x-x_0)\delta(y-y_0) \sin\theta\,\mathrm{d}\theta \,\mathrm{d}\phi. $$
I know that this integral is very simple when the delta function is fully three-dimensional, such as $$\delta^3(\mathbf{r}-\mathbf{r}_0)=\delta(x-x_0)\delta(y-y_0)\delta(z-z_0),$$ but what I have here is something a bit different. Does anyone have anyone have any suggestions?
Here we give a rough sketch, and leave it to the reader to iron out various technical issues.
Define the spherical coordinates $(r,\theta,\phi)$ as $$\begin{align} x~=~&r\xi, \qquad \xi~=~\sin\theta\cos\phi, \cr y~=~&r\eta, \qquad \eta~=~\sin\theta\sin\phi, \cr z~=~&r\cos\theta. \end{align} \tag{1} $$
Define reflection $$\sigma: (x,y,z)\mapsto (x,y,-z)\tag{2}$$ in the $(x,y)$-plane.
Define function $$g({\bf r})~:=~f({\bf r})+f\circ\sigma({\bf r}), \qquad {\bf r}~=~(x,y,z).\tag{3}$$
Define function $$h({\bf r})~:=~\frac{g({\bf r})}{rz}.\tag{4}$$
Define unit disk $$ D~:=~\{(\xi,\eta)\in\mathbb{R}^2 \mid \xi^2+\eta^2~<~1\}. \tag{5}$$
OP's integral then reads $$\begin{align} I:~=~& \int_0^{2\pi}\!d\phi~\int_0^{\pi}\!d\theta~\sin\theta ~f({\bf r}) \delta(x-x_0)\delta(y-y_0)\cr ~=~&\int_0^{2\pi}\!d\phi~\int_0^{\frac{\pi}{2}}\!d\theta~\sin\theta\cos\theta ~h({\bf r})\delta\left(\xi-\frac{x_0}{r}\right)\delta\left(\eta-\frac{x_0}{r}\right)\cr ~=~&\iint_D \!d\xi~d\eta ~h({\bf r})\delta\left(\xi-\frac{x_0}{r}\right)\delta\left(\eta-\frac{x_0}{r}\right)\cr ~=~&\left\{\begin{array}{ccc} h\left(x_0, y_0, \sqrt{r^2-x_0^2-y_0^2}\right) &{\rm for}& x_0^2+y_0^2<r^2, \cr \frac{1}{2}h(x_0, y_0, 0) &{\rm for}& x_0^2+y_0^2=r^2, \cr 0 &{\rm for}& x_0^2+y_0^2>r^2. \end{array}\right. \end{align}\tag{6}$$