Evaluating $\int_0^L\int_0^L \cos^{-1}\left(\frac{x^2+y^2-a^2}{2xy}\right)xy~dx~dy$ I have the following integral to evaluate.

$$\displaystyle\int\int_{x,y\in\Omega}\cos^{-1}\left(\frac{x^2+y^2-a^2}{2xy}\right)xy~dx~dy,$$
  where $0\le a\le L$ and $\Omega=\{(x,y)\in (0,L]\times(0,L]: |x-y|<a, 
(x+y)>a\}$.

I used the transformation of variable as below: $$\begin{align}x&=\frac{a\sin(\theta+\phi)}{\sin \theta},\\
\\
y&=\frac{a\sin\phi}{\sin \theta}.\end{align}$$ It seems that the integral can be evaluated by this technique but I'm not sure because the calculations are getting quite messy. I was wondering is there any better method to find it? Thanks in advance. 
 A: The integration domain is separated in four parts (figure below). The analytical result is obtained for parts (1) and (2) where the calculus is more difficult : result fonction of two paraleters $(a, L)$. The parts 3 and 4 are simpler : one parameter only $(a)$. It is solved with the same method.

A: I'd much rather keep the integration domain square.
My suggestion is to take the derivative over parameter $a$. This makes the integral much more manageable. Arccos is gone, the square root simplifies and the $xy$ vanishes. That integral should be easier to evaluate. Then, integrate in $a$ to retrieve the original answer.
The initial value at $a=0$ still has to be integrated by hand though.
EDIT:
$$F(a)=\iint \cos^{-1}\left(\frac{x^2+y^2-a^2}{2xy}\right)xy{\,\rm d}x{\,\rm d}y$$
Derivative with respect to $a$:
$$F'(a)=\iint \frac{-1}{\sqrt{1-\left(\frac{x^2+y^2-a^2}{2xy}\right)^2}}\left(\frac{-a}{xy}\right)xy{\,\rm d}x{\,\rm d}y$$
$$F'(a)=a\iint \frac{1}{\sqrt{\frac{4x^2y^2-\left(x^2+y^2-a^2\right)^2}{(2xy)^2}}}{\,\rm d}x{\,\rm d}y$$
$$F'(a)=a\iint \frac{2xy}{\sqrt{2a^2(x^2+y^2)-a^4-(x^2-y^2)^2}}{\,\rm d}x{\,\rm d}y$$
This may be easier to integrate. At least once, it can be done almost trivially ($u=x^2$). Or, you could use $u=x^2+y^2$ and $v=x^2-y^2$, which sets the integration domain to the diagonal square bound by $u+v=2L$, $u-v=2L$, $u+v=0$, $u-v=0$ which you could cut in two for symmetry. That's not much better, but at least less writing.
$$F'(a)=a\left(\int_0^{L}\int_0^{u}+\int_L^{2L}\int_0^{2L-u}\right) \frac{1}{\sqrt{2a^2u-a^4-v^2}}{\,\rm d}v{\,\rm d}u$$
Note that this is just a suggestion - the second integation is still horrible and may not yield results. I don't have Mathematica with me so I'm not checking now.
If you manage to evaluate the integral now, you just integrate again:
$$F(a)=F(0)+\int_0^a F(t){\,\rm d}t$$
