integration - How to calculate this double integral? (solution known)

This question is related to determining the one-sided spectral density function $G(\omega)$ of a 2D random field which is given by:

$$G(\omega_1,\omega_2) = \left(\frac{2 \sigma}{\pi}\right)^2 \int_{0}^{\infty} \int_{0}^{\infty} \exp\left(-\sqrt{\left(\frac{2\tau_1}{\theta_1}\right)^2 + \left(\frac{2\tau_2}{\theta_2}\right)^2}\right) \cos(\omega_1 \tau_1) \cos(\omega_2 \tau_2)\ \, d\tau_1 d\tau_2$$

Fenton and Griffiths (2008) in 'Risk Assessment in Geotechnical Engineering' give the analytical result as:
\begin{equation}
G(\omega_1,\omega_2) = \frac{\sigma^2 \theta_1 \theta_2}{2 \pi \left[1 + \left(\frac{\theta_1 \omega_1}{2}\right)^2 + \left(\frac{\theta_2 \omega_2}{2}\right)^2 \right]^{3/2}}
\end{equation}

However I have been unable to obtain this result myself. I have also tried integration in Wolfram Mathematica which failed as well.

Does anyone know how to perform this integral or formulate it in such a way that Wolfram Mathematica can solve it? I want to verify the solution given by Fenton and Griffiths (2008).

Answer

Hint. One may start with the classic result
$$
\int_0^\infty e^{- \lambda x}\:dx=\frac1{\lambda},\qquad \text{Re}(\lambda) >0, \tag1
$$ giving with $\lambda:= 1+ia$ and by differentiating
$$

\int_0^\infty e^{-r}\cdot r \cdot \cos(ar)\:dr=\frac{1-a^2}{\left(1+a^2\right)^2},\qquad a \in \mathbb{R}. \tag2
$$ Then, by a standard change of variable from cartesian coordinates to polar coordinates, one gets
$$
\begin{align}
&G(\omega_1,\omega_2) = \left(\frac{2}{\pi}\right)^2\!\! \int_{0}^{\infty} \!\!\int_{0}^{\infty}\!\! \exp\left(-\sqrt{\left(\frac{2\tau_1}{\theta_1}\right)^2 + \left(\frac{2\tau_2}{\theta_2}\right)^2}\right) \cos(\omega_1 \tau_1) \cos(\omega_2 \tau_2) \: d\tau_1 d\tau_2
\\&=\frac{\theta_1\theta_2}{\pi^2}\!\! \int_{0}^{\infty} \!\!\int_{0}^{\infty}\!\! \exp\left(-\sqrt{x^2 + y^2}\right) \cos(w_1 x) \cos(w_2 y) \: dx dy
\\&=\frac{\theta_1\theta_2}{\pi^2}\!\! \int_{0}^{\pi/2} \!\!\int_{0}^{\infty}\!\! re^{-r}\cos(w_1 r \cos \theta) \cos(w_2 r \sin \theta) \: dr d\theta
\\&=\frac{\theta_1\theta_2}{2\pi^2}\!\! \int_{0}^{\pi/2} \!\!\int_{0}^{\infty}\!\! re^{-r} \left(\cos[(w_1 \cos \theta+w_2 \sin \theta)r]+ \cos[(w_1 \cos \theta-w_2 \sin \theta)r]\right) dr d\theta
\\&=\frac{\theta_1\theta_2}{2\pi^2}\!\! \int_{0}^{\pi/2} \!\!\left(\frac{1-(w_1 \cos \theta+w_2 \sin \theta)^2}{\left(1+(w_1 \cos \theta+w_2 \sin \theta)^2\right)^2}+\frac{1-(w_1 \cos \theta-w_2 \sin \theta)^2}{\left(1+(w_1 \cos \theta-w_2 \sin \theta)^2\right)^2}\right) d\theta
\\&=\frac{\theta_1\theta_2}{2\pi^2}\cdot \frac{\pi}{(1+w_1^2+w_2^2)^{3/2}}

\end{align}
$$ where we have obtained the latter integral by a tangent half-angle substitution and where we have set
$$
w_1:=\frac{\theta_1 \omega_1}2,\quad w_1:=\frac{\theta_2 \omega_2}2.
$$