Solution of the integral $\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\frac{3xy}{(x^{2}+y^{2}+z^{2})^{5/2})}e^{i(k_{x}x+k_{y}y)} dx dy$ I'm trying to solve the following integral:
$$\int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty}\frac{3xy}{(x^{2}+y^{2}+z^{2})^{5/2})}e^{i(k_{x}x+k_{y}y)} dx dy$$
i'm using the solution posted as the answer to another question: here's the link
The principle behind the solution is to change the coordinates of integration and obtain a separation of the integrals.
In particular following the steps reported in the link. The exponential argument can be written as $e^{\boldsymbol{k}\cdot \boldsymbol{x}}$, where $\boldsymbol{k}= k_{x} \hat{x} +k_{y} \hat{y}$ moreover $\boldsymbol{k} = k \hat{k}$ here we have defined $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$. 
As reported in the link we can orient $\hat{k}$ such that  $\mathbf{k}\mapsto k\hat{x}$.
Considering also $s=\sqrt{x^{2}+y^{2}}$ and substituting the new variables in the integral I obtain:
$\int_{0}^{\infty}\int_{-\pi}^{\pi}\frac{Cos(\theta)Sin(\theta) \; e^{ik s\cos\theta }}{\sqrt{s^2+z^2}}s^3\,d\theta \,ds$
that can be integrated with mathematica for example returning zero because the integral:
$\int_{-\pi}^{\pi}Cos(\theta)Sin(\theta) \; e^{ik s\cos\theta }d\theta=0$
however the paper that i'm following reports the following as solution of the integral:
$-\frac{k_{x}k_{y}}{k}e^{-kz}$
 A: Here is an answer in the spirit of my second approach of your linked question. We want to find the planar Fourier transform of $F(\mathbf{x})=\dfrac{3xy}{(x^2+y^2+z^2)^{5/2}}$. Let's start with the 3D Fourier transform. In the previous problem, we were able to take advantage of the fact that the function $G(\mathbf{x})=(x^2+y^2+z^2)^{-1/2}$ was the Green's function of the Laplacian i.e. $\nabla^2 G=-4\pi \delta^{3}(\mathbf{r})$ and thus had the simple Fourier transform $(k_x^2+k_y^2+k_z^2)G(\mathbf{k})=4\pi$.
Happily, we can take advantage of this again: observe that $$\frac{\partial^2 G}{\partial y\partial x} =\frac{\partial}{\partial x}\frac{-y}{(x^2+y^2+z^2)^{3/2}}=\frac{3xy}{(x^2+y^2+z^2)^{5/2}}=F(\mathbf{x}),$$ so the Fourier transform simply gives
$$ F(\mathbf{k})=-k_x k_y G(\mathbf{k})=- \dfrac{4\pi k_x k_y}{k_x^2+k_y^2+k_z^2}.$$
Even better, this only differs from my last answer in the $k_x,k_y$-dependence; the $z$-dependence is unchanged, and so the same result for the inverse Fourier transform in $z$ is applicable. Consequently 
$$\boxed{\displaystyle F(k_x,k_y,z)=-2\pi k_x k_y \left(\dfrac{e^{-|z|\sqrt{k_x^2+k_y^2}}}{\sqrt{k_x^2+k_y^2}}\right)},$$ which apart from a factor of $2\pi$ due to definitions is exactly the expected result.
