I am trying to analytically evaluate
$$f(a,b)=\int_0^{2\pi} K_0(\sqrt{a^2+b^2-2ab\cos(\theta)}) d \theta$$
where $K_0$ is a modified Bessel function of the second kind and $a,b>0$. I happen to know the solution by a different method, it is
$$f(a,b)=2\pi I_0(\min \{ a,b \}) K_0(\max \{ a,b \})$$
where $I_0$ is a modified Bessel function of the first kind. But I don't know how to evaluate the integral directly, and CAS's seem to choke on it (though they can check the result above numerically just fine). It seems like it is probably useful to use the identity
$$K_0(x)=\int_0^\infty \frac{\cos(x t)}{\sqrt{t^2+1}} dt$$
and then interchange the order of integration, but then I still need to be able to compute
$$\int_0^{2\pi} \cos(\sqrt{a^2+b^2-2ab\cos(\theta)}) d \theta$$
for arbitrary $a,b>0$, before proceeding to the $dt$ integral.
Note that the problem is significantly simpler when $a=b$ because in this case $\sqrt{a^2+b^2-2ab\cos(\theta)}=2a\sin(\theta/2)$; in this case the cosine integral is $2\pi J_0(2at)$, which results in the easier problem
$$2\pi \int_0^\infty \frac{J_0(2at)}{\sqrt{t^2+1}}.$$
But in the case $a \neq b$ I don't even see how to begin.
This is a consequence of the rather miraculous-looking Graf addition formula: if we assume that $0<b<a$ (since the symmetry is apparent), and write $w=\sqrt{a^2+b^2-2ab\cos{\theta}}$, we have $$ J_{\nu}(w) = \left(\frac{a-be^{i\theta}}{a-be^{-i\theta}}\right)^{\nu/2} \sum_{m=-\infty}^{\infty} J_{\nu+m}(a)J_{m}(b) e^{im\theta}. $$ One can prove this directly using the integral representation of $J_{\nu}$ on the right several times and swapping various integrals and sums: the proof is given in Watson's book, p. 360). Mucking about with this in the standard way gives the corresponding formulae for modified Bessel functions, (also given in Watson's book, on the next page) of which the one of interest to us is $$ K_{0}(w) = \sum_{m=-\infty}^{\infty} K_m(a) I_m(b) e^{im\theta}. $$ Integrating this from $0$ to $2\pi$ gives the result since only the $m=0$ term survives.