Integration of $\int_0^\pi\int_0^\pi[\sin x\sin y\ge a]\,dx\,dy$ Does the integral $$I(a)=\int_0^\pi\int_0^\pi[\sin x\sin y\ge a]\,dx\,dy$$ have any closed-form solution? This is the area under the contours of $\sin x\sin y$; $[\cdot]$ is the Iverson bracket, which is the indicator function for the predicate in the bracket. (For the curious, this integral was the result of a thought experiment about how to make those dots used for shading in comics to produce a given brightness level.)
Obviously, $I(0)=\pi^2$ and $I(1)=0$, and $I'(1)=-2\pi$ because if $a=1-\epsilon$, the approximation $\sin x\approx 1-\frac12(x-\frac\pi2)^2$ suggests that the area is well approximated by a circle $\frac12(x-\frac\pi2)^2+\frac12(y-\frac\pi2)^2\le\epsilon$, with area $2\pi\epsilon$. Are there any other conclusions that can be drawn? What is $I'(0)$?
EDIT: $I'(0)$ is $-\infty$, and the best low-order approximation for $I$ near $0$ is $I(a)\approx\pi^2-4a(\log\frac4a+1)$. For $a=\epsilon$, the area "missing" from the box $\pi^2$ can be divided into the four corners (the largest squares that can fit into the four corners of the box, of side length $\sin^{-1}\sqrt a$ and area $\approx a$ each), plus the four sides, which are given by the integral $$8\int_\sqrt a^{\pi/2}\frac a{\sin x}\,dx=8a\log\cot\frac{\sqrt a}2\approx8a\log\frac2{\sqrt a}=4a\log\frac4a,$$ hence the approximation $I(a)\approx\pi^2-4a(\log\frac4a+1)$.
 A: By symmetry, $I(a)$ is $4$ times the same integral but on $(0,\frac\pi2)^2$ only. The change of variable $$(z,s)=(\sin x\sin y,\cos x),$$ sends the domain $(0,\frac\pi2)^2$ to the positive quarter of the unit disk defined by $$z\gt0,\quad s\gt0,\quad z^2+s^2\lt1.$$ The Jacobian is given by $\mathrm dz\mathrm ds=\sin^2x\cos y\mathrm dx\mathrm dy$ and $$\sin^2x\cos y=\sqrt{1-s^2}\cdot\sqrt{1-z^2-s^2},$$ hence
$$
I(a)=4\int_a^1\mathrm dz\int_0^{\sqrt{1-z^2}}\frac{\mathrm ds}{\sqrt{1-s^2}\cdot\sqrt{1-z^2-s^2}}.
$$
The change of variable $s=\sqrt{1-z^2}\cdot t$ yields finally
$$
I(a)=4\int_a^1\mathrm dz\int_0^1\frac{\mathrm dt}{\sqrt{1-t^2}\cdot\sqrt{1-(1-z^2)t^2}},
$$
that is,
$$
I(a)=4\int_a^1K(\sqrt{1-z^2})\mathrm dz,
$$
where $K$ denotes the complete elliptic integral of the first kind. Thus, for every $a$ in $(0,1]$,
$$
I'(a)=-4K(\sqrt{1-a^2}),
$$
for example $I'(1)=-4K(0)=-2\pi$ and, at $a=0^+$, $
I'(0^+)=-\infty.
$ 
An equivalent formula, based on the hypergeometric function ${}_2F_1$, is
$$
I(a)=2\pi\int_a^1{}_2F_1(\tfrac12,\tfrac12;1;1-z^2)\mathrm dz,
$$
