How do I calculate the triple integral $\int_{0}^{2\pi}\int_0^1 \int_0^1 xy\sqrt{x^2 + y^2 -2xy\cos(\theta)} dx \text{ }dy \text{ } d\theta$? I have tried a couple substitutions which didn't pan out, and I tried differentiating under the integral sign on
$$
I(\theta) = \int_0^1 \int_0^1 xy\sqrt{x^2 + y^2 -2xy\cos(\theta)} dx \text{ }dy
$$
to get
$$
\frac{dI}{d\theta} = -\sin(\theta) \int_0^1 \int_0^1 \frac{x^2y^2}{\sqrt{x^2 + y^2 -2xy\cos(\theta)}} dx \text{ }dy
$$
but nothing has worked so far. Can any of you help me?
 A: By symmetry, we have:
$$J=\int_{0}^{2\pi}I(\theta)\,d\theta = 2\int_{0}^{2\pi}\iint_{D}xy\sqrt{x-ye^{i\theta}}\sqrt{x-ye^{-i\theta}}\,dx\,dy\,d\theta $$
where $D$ is the region:
$$ D = \{(x,y)\in[0,1]^2 : y\leq x\},$$
so:
$$\begin{eqnarray*} J &=& 2\int_{0}^{2\pi}\int_{0}^{1}\int_{0}^{x}xy\sqrt{x-ye^{i\theta}}\sqrt{x-ye^{-i\theta}}\,dy\,dx\,d\theta\\&=&2\int_{0}^{2\pi}\int_{0}^{1}x^4\int_{0}^{1}t\sqrt{1-te^{i\theta}}\sqrt{1-te^{-i\theta}}\,dt\,dx\,d\theta\\&=&\frac{2}{5}\int_{0}^{1}t\int_{0}^{2\pi}\sqrt{1-te^{i\theta}}\sqrt{1-te^{-i\theta}}\,d\theta\,dt.\end{eqnarray*}$$
Since for any $|z|<1$ we have:
$$\sqrt{1-z}=\sum_{n\geq 0}\binom{\frac{1}{2}}{n}(-1)^n z^n $$
it happens that:
$$ \int_{0}^{2\pi}\sqrt{1-te^{i\theta}}\sqrt{1-te^{-i\theta}}\,d\theta =2\pi\sum_{n\geq 0}\binom{\frac{1}{2}}{n}^2 t^{2n} $$
so:
$$ J = \frac{4\pi}{5}\int_{0}^{1}\sum_{n\geq 0}\binom{\frac{1}{2}}{n}^2 t^{2n+1}\,dt = \frac{4\pi}{5}\sum_{n\geq 0}\binom{\frac{1}{2}}{n}^2\frac{1}{2n+2}=\frac{4\pi}{5}\cdot\frac{16}{9\pi}=\color{red}{\frac{64}{45}}.$$
To prove the last identity, notice that:
$$\sum_{n\geq 0}\binom{\frac{1}{2}}{n}^2\frac{1}{2n+2}=\sum_{n\geq 0}\frac{\frac{1}{16^n}\binom{2n}{n}^2}{(2n-1)^2(2n+2)}$$
can be computed from:
$$ \sum_{n\geq 0}\frac{\frac{1}{4^n}\binom{2n}{n}}{(2n-1)^2}x^{2n}=\sqrt{1-x^2}+x\arcsin x,$$
$$ \sum_{n\geq 0}\frac{\frac{1}{4^n}\binom{2n}{n}}{(2n+2)}x^{2n}=\frac{1-\sqrt{1-x^2}}{x^2}$$
through a contour integral, as already guessed by the OP.

This can be tackled also in an alternative and very elegant way. 
If we consider the extremely special function:
$$ B(\lambda)=\sum_{n\geq 0}\left(\frac{\binom{2n}{n}}{(2n-1)4^n}\right)^2\lambda^{2n}=\frac{1}{2\pi}\int_{0}^{2\pi}\sqrt{1+\lambda^2-2\lambda\cos\theta}\,d\theta=\phantom{}_2 F_1\left(-\frac{1}{2},-\frac{1}{2};1;\lambda^2\right)$$
it is easy to check from the power series representation that $B$ satisfies the ODE:
$$ \lambda B = (\lambda^2+1) B' + (\lambda-\lambda^3) B'' $$
but from that differential equation it is also easy to check, by integration by parts, that:
$$ \int_{0}^{1} \lambda\, B(\lambda)\,d\lambda = \frac{16}{9\pi} $$
since $B(1)=\frac{4}{\pi}=2\cdot B'(1)$. The previous integral is everything we need to solve our problem.
