Can the following triple integral be computed via elementary calculus methods?

Consider the following triple integral:

$$\int_0^{2\pi}\int_0^1 \int_0^1 xy\sqrt{x^2 + y^2 -2xy\cos(\theta)} \, dx \, dy \, d\theta$$

A solution was provided to this integral by Jack D'Aurizio here, but both his solutions required rather sophisticated methods, like elliptic integrals and special function expansions. What I've been struggling with-and would like a second opinion on-is whether or not this integral can be solved in closed form using very simple calculus techniques,like a standard change of variables to plane polar in the xy domain or spherical coordinates in $R^3$.

My labor over the last 2 days,multiple false starts and geometric arguments in the domain seem to indicate that the answer is no because there's no way to set up the integral without introducing a term of $\sqrt {\sin (ax)}$ or $\sqrt {\cos(ax)}$ at some point. Therefore, some special functional substitution or numerical method solution is needed.

Or am I wrong?

Answer

$$
\begin{align}
&\int_0^{2\pi}\int_0^1\int_0^1xy\sqrt{x^2+y^2-2xy\cos(\theta)}\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}\theta\tag{1}\\
&=2\int_0^{2\pi}\int_0^1\int_0^yxy\sqrt{x^2+y^2-2xy\cos(\theta)}\,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}\theta\tag{2}\\
&=2\int_0^{2\pi}\int_0^1\int_0^1ry^2\sqrt{r^2y^2+y^2-2ry^2\cos(\theta)}\,\,y\,\mathrm{d}r\,\mathrm{d}y\,\mathrm{d}\theta\tag{3}\\
&=2\int_0^{2\pi}\int_0^1\int_0^1ry^4\sqrt{r^2+1-2r\cos(\theta)}\,\mathrm{d}y\,\mathrm{d}r\,\mathrm{d}\theta\tag{4}\\
&=\frac25\int_0^{2\pi}\int_0^1\sqrt{r^2+1-2r\cos(\theta)}\,r\,\mathrm{d}r\,\mathrm{d}\theta\tag{5}\\

&=\frac25\int_{-\pi/2}^{\pi/2}\int_0^{2\cos(\theta)}r^2\,\mathrm{d}r\,\mathrm{d}\theta\tag{6}\\
&=\frac2{15}\int_{-\pi/2}^{\pi/2}8\cos^3(\theta)\,\mathrm{d}\theta\tag{7}\\
&=\frac{16}{15}\int_{-1}^1(1-u^2)\,\mathrm{d}u\tag{8}\\[4pt]
&=\frac{64}{45}\tag{9}
\end{align}
$$
Explanation:
$(2)$: the integral is the same for $x\lt y$ as for $x\gt y$, so assume $x\lt y$ and multiply by $2$
$(3)$: substitute $r=\frac xy$
$(4)$: collect the $y$s and switch the order of integration
$(5)$: integrate in $y$
$(6)$: $(5)$ is the distance from $(1,0)$ integrated over the unit disk centered at $(0,0)$;
$\hphantom{(6):}$this is the same as $r$ integrated over the unit disk centered at $(1,0)$,
$\hphantom{(6):}$whose equation is $r\le2\cos(\theta)$ for $\theta\in[-\pi/2,\pi/2]$
$(7)$: integrate in $r$
$(8)$: substitute $u=\sin(\theta)$
$(9)$: integrate in $u$