Evaluating $\int_0^1\int_0^1\left(\sqrt{1-x^2}+\sqrt{1-y^2}\right)^{-1}\,dx\,dy$ I am trying to get $$\mathcal{J}=\int_0^1\int_0^1\frac{1}{\sqrt{1-x^2}+\sqrt{1-y^2}} \,dx\,dy,\tag{1}$$
as variation of the integral $\int_0^1\frac{dt}{\sqrt{1-t^2}}$. Then exploiting the same change of varible, say us $x=\cos u$ that we use in this last integral we get $$\mathcal{J}=\int_0^1dy\int_0^{\pi/2}\frac{\sin u}{\sqrt{1-y^2}+\sin u}\,du\tag{2}.$$
We denote $A=\sqrt{1-y^2}$, and thus $\sqrt{A^2-1}=-y^2$. Then I know uisng a CAS that the integral 
$$\int_0^{\pi/2}\frac{\sin u}{A+\sin u}\,du$$
equals to $$\frac{\pi}{2}-\frac{2A\tan^{-1}\left(\frac{1+A}{\sqrt{A^2-1}}\right)}{\sqrt{A^2-1}}+\frac{2A\tan^{-1}\left(\frac{1}{\sqrt{A^2-1}}\right)}{\sqrt{A^2-1}}.$$

Question. How can we finish the calculations (if it is possible) to get the closed-form of $$\mathcal{J}=\int_0^1\int_0^1\frac{1}{\sqrt{1-x^2}+\sqrt{1-y^2}}\,dx\,dy?$$
  Many thanks.

I'm interested in some evaluation of such integral, thus if it isn't possible to get a closed-form show your approach.
 A: Let's see. We have
$$\mathcal{J}=\iint_{(0,\pi/2)^2}\frac{\cos\theta\cos\varphi}{\cos\theta+\cos\varphi}\,d\theta\,d\varphi=\frac{1}{2}\int_{0}^{\pi/2}\cot\theta\left(\pi\sin\theta-4\cos\theta\,\text{arctanh}\tan\frac{\theta}{2}\right)\,d\theta$$
and the RHS is related to the inverse Gudermannian function.
$$\begin{eqnarray*}\mathcal{J}&=&\int_{0}^{\pi/4}\cot(2\theta)\left(\pi\sin(2\theta)-4\cos\theta\,\text{arctanh}\tan\theta\right)\,d\theta\\&=&\int_{0}^{1}\frac{1-x^2}{2x(1+x^2)}\left(\frac{2\pi x}{1+x^2}-\frac{4(1-x^2)\text{arctanh}(x)}{1+x^2}\right)\,dx\\&=&\frac{\pi}{2}-2\int_{0}^{1}\left(\frac{1-x^2}{1+x^2}\right)^2\frac{\text{arctanh}(x)}{x}\,dx\end{eqnarray*}$$
and by integration by parts the problem boils down to evaluating $\int_{0}^{1}\frac{\log(x)}{1-x^2}\,dx=-\frac{\pi^2}{8}$:
$$\boxed{\mathcal{J}= \frac{\pi(4-\pi)}{4}.}$$
$\mathcal{J}$ can also be represented as
$$ \int_{0}^{+\infty}\left(\int_{0}^{\pi/2}\cos\theta e^{-x\cos\theta}\,d\theta\right)^2\,dx =\frac{\pi^2}{4}\int_{0}^{+\infty}\left(L_{-1}(x)-I_1(x)\right)^2\,dx$$
in terms of a Struve function and a modified Bessel function of the first kind. $L_{-1}(x)-I_1(x)$ is very regular on $\mathbb{R}^+$: it approximately behaves like $\frac{2}{\pi} e^{-\pi x/4}$, immediately leading to $\mathcal{J}\approx \frac{2}{\pi}$. Numerically
$$ \mathcal{J} = 0.674191\ldots \approx \frac{2001}{2968}. $$

Folklore: a nice consequence of $\frac{\pi(4-\pi)}{4}<\frac{2}{\pi}$ is $\color{blue}{\pi<1+\sqrt{5}=2\varphi}$.
A: The exact answer is 
$$ \mathcal{J}=\int_0^1 \int_0^1 \frac{dx \, dy} {\sqrt{1-x^2}+\sqrt{1-y^2}} =
\frac{\pi}{4}(4-\pi).$$
I derived it by letting a CAS do the alternate integral for me:
$$ \mathcal{J} =
\frac{2}{\pi}\int_0^\infty\Big( \frac{u}{\sqrt{1+u^2}}\,\tanh^{-1}(1/\sqrt{1+u^2} \Big)^2.$$
A sketch of a proof of the ID is as follows.  Use the Euler gamma representation to show 
$$ \mathcal{J} = \int_0^\infty dt \Big( \int_0^1 dx\,\exp{(-t\sqrt{1-x^2})}
\Big)^2. $$
For each inner integral use the ID
$$ \int_0^\infty \frac{du \,u}{b^2+u^2}\sin{(a \,u)} = \frac{\pi}{2}
 \exp{(-a\, b)} $$ so you have a four-fold integral.
Switch integral signs.  You end up with an inner integral that looks like
$$ \int_0^\infty \sin{(t \,u)} \sin{(t \,w)} dt,$$
which is related to the Dirac delta function.  4-fold integral collapses to a 3-fold, and the two innermost ones separate and have the exact form as given in the integrand of the second formula above.
