Month ago I encounter a nice result numerically checked by Mathematica $$ \int_{0}^{\pi/2} \int_{0}^{\pi/2} \frac{\theta\cot\theta-\varphi\cot\varphi}{\cos\theta-\cos\varphi} \mathrm{d}\varphi\mathrm{d}\theta = \pi\ln2 $$ where the integrated function is actually well-defined even around its singularity $\theta=\varphi=0$.
At my first sight, I thought it might be a trivial conclusion derived from a kind of typical integral like $$ \int_{0}^{\pi} \frac{\cos n\theta}{\cos\theta-\cos\varphi} \mathrm{d}\theta = \pi\frac{\sin n\varphi}{\sin\varphi} $$ just using a proper series expansion. However, when I review it in detail, the result over $(0,\pi/2)$ will be awkwardly complicated. I realize this double integrals may not be done directly, or I may lack some essential insight to solve it.
So I question it here for some further discussion, and thanks in advance for any suggestion.

Here is a solution based on @G Cab's remarkable observation.
Reduction. Let $I$ denote the integral, and apply the substitution $(\theta,\varphi)\mapsto(2\theta,2\varphi)$ to write
$$ I = 4 \int_{0}^{\pi/4} \int_{0}^{\pi/4} \frac{2\theta\cot(2\theta) - 2\varphi \cot(2\varphi)}{\cos(2\theta)-\cos(2\varphi)} \, \mathrm{d}\theta\mathrm{d}\varphi. $$
If we define the $1$-form $\omega$ by
$$ \omega = \left( \log\frac{\cot\theta}{\cot\varphi} \right) \biggl( \frac{\mathrm{d}( (\theta-\varphi)^2 )}{2\sin^2(\theta-\varphi)} - \frac{\mathrm{d}( (\theta+\varphi)^2 )}{2\sin^2(\theta+\varphi)} \biggr), $$
then $\omega$ is smooth in the region $(0, \pi/2)^2$ via continuation along the diagonal $\varphi = \theta$. Moreover, @G Cab's computation for the differential reduces to1)
$$ \mathrm{d}\omega = 4 \left( \frac{2\theta\cot(2\theta) - 2\varphi \cot(2\varphi)}{\cos(2\theta)-\cos(2\varphi)} \right) \, \mathrm{d}\theta \wedge \mathrm{d}\varphi. $$
Now let $\epsilon > 0$ be small and consider the triangular region
$$ \mathcal{T}_{\epsilon} = \bigl\{ (\theta, \varphi) : \epsilon \leq \varphi \leq \theta \leq \tfrac{\pi}{4} \bigr\}. $$
Its boundary $\partial\mathcal{T}_{\epsilon}$ consists of three oriented line segments
Then by the symmetry and the Stokes' theorem,
\begin{align*} I &= 2 \lim_{\epsilon \to 0^+} \int_{\mathcal{T}_{\epsilon}} \mathrm{d}\omega = 2 \lim_{\epsilon \to 0^+} \int_{\partial \mathcal{T}_{\epsilon}} \omega = 2 (I_{\rightarrow} + I_{\uparrow} + I_{\swarrow}), \end{align*}
where $\rightarrow$, $\uparrow$, and $\swarrow$ denote the oriented bottom/height/hypothenus of $\partial\mathcal{T}_{\epsilon}$, respectively, and
$$ I_{\bullet} = \lim_{\epsilon \to 0^+} \int_{\bullet} \omega $$
for each $\bullet \in \{ \rightarrow, \uparrow, \swarrow \}$.
Integral Computations. We will evaluate each $I_{\bullet}$:
1. It is not hard to check that $\omega = 0$ along $\varphi = \theta$, and so, we get $I_{\swarrow} = 0$.
2. For $I_{\uparrow}$, we note that
\begin{align*} I_{\uparrow} &= \int_{\varphi = 0}^{\varphi = \pi/4} \omega|_{\theta=\pi/4} \\ &= \int_{\varphi = 0}^{\varphi = \pi/4} \left( - \log\cot\varphi \right) \biggl( \frac{\mathrm{d}( (\pi/4-\varphi)^2 )}{2\sin^2(\pi/4-\varphi)} - \frac{\mathrm{d}( (\pi/4+\varphi)^2 )}{2\sin^2(\pi/4+\varphi)} \biggr) \\ &= 2 \int_{0}^{\pi/4} \left( \log\cot\varphi \right) \biggl( \frac{\pi/4-\varphi}{1-\sin(2\varphi)} + \frac{\pi/4+\varphi}{1+\sin(2\varphi)} \biggr) \, \mathrm{d}\varphi \\ &= 2 \int_{-\pi/4}^{\pi/4} \frac{\pi/4+\varphi}{1+\sin(2\varphi)} \log \left| \cot\varphi \right| \, \mathrm{d}\varphi \\ &= \frac{1}{2} \int_{0}^{\pi} \frac{x}{1-\cos x} \log \left| \cot \left(\frac{x}{2} - \frac{\pi}{4}\right) \right| \, \mathrm{d}x, \end{align*}
where the last step is the result of the substitution $\varphi=\tfrac{x}{2}-\tfrac{\pi}{4}$. Then by noting that
$$ \log \left| \cot \left(\frac{\pi}{4} - \frac{x}{2} \right) \right| = \operatorname{Re}\log\left(\frac{1-ie^{ix}}{1+ie^{ix}}\right) = i (\arctan(e^{-ix}) - \arctan(e^{ix})), $$
we get
\begin{align*} I_{\uparrow} &= \frac{1}{2}\int_{0}^{\pi} \frac{ix}{1-\cos x} (\arctan(e^{-ix}) - \arctan(e^{ix})) \, \mathrm{d}x \\ &= \frac{1}{2}\int_{-\pi}^{\pi} \frac{ix}{1-\cos x} \left(\frac{\pi}{4} - \arctan(e^{ix})\right) \, \mathrm{d}x \\ &= i \int_{\mathcal{C}} \frac{(\log z )(\frac{\pi}{4} - \arctan z)}{(z-1)^2} \, \mathrm{d}z, \end{align*}
where the substitution $z = e^{ix}$ is applied to the last step and $\mathcal{C}$ is the path parametrized by $e^{ix}$ for $-\pi < x < \pi$. Now deforming the contour $\mathcal{C}$ so as to enclose the principal branch cut for $\log$,
\begin{align*} I_{\uparrow} &= 2\pi \int_{0}^{1} \frac{\frac{\pi}{4} + \arctan t}{(t+1)^2} \, \mathrm{d}t = 2\pi \int_{0}^{1} \frac{1}{(t^2+1)(t+1)} \, \mathrm{d}t = \frac{\pi^2}{4} + \frac{\pi}{2} \log 2. \end{align*}
3. For the integral along the bottom,
\begin{align*} \int_{\rightarrow} \omega &= \int_{\theta = \epsilon}^{\theta = \pi/4} \omega|_{\varphi=\epsilon} \\ &= \int_{\theta = \epsilon}^{\theta = \pi/4} \log\left( \frac{\cot\theta}{\cot\epsilon} \right) \biggl( \frac{\mathrm{d}( (\theta-\epsilon)^2 )}{2\sin^2(\theta-\epsilon)} - \frac{\mathrm{d}( (\theta+\epsilon)^2 )}{2\sin^2(\theta+\epsilon)} \biggr) \\ &= \int_{\epsilon}^{\pi/4} \log\left( \frac{\cot\theta}{\cot\epsilon} \right) \biggl( \frac{\theta-\epsilon}{\sin^2(\theta-\epsilon)} - \frac{\theta+\epsilon}{2\sin^2(\theta+\epsilon)} \biggr) \, \mathrm{d}\theta. \end{align*}
Substituting $\theta = \epsilon x$,
\begin{align*} \int_{\rightarrow} \omega &= \int_{1}^{\pi/4\epsilon} \log\left( \frac{\cot(\epsilon x)}{\cot\epsilon} \right) \biggl( \frac{\epsilon^2 (x-1)}{\sin^2(\epsilon(x-1))} - \frac{\epsilon^2 (x+1)}{2\sin^2(\epsilon(x+1))} \biggr) \, \mathrm{d}x. \end{align*}
It is not hard (although not entirely trivial) to show that the above integral converges to
\begin{align*} I_{\rightarrow} = \lim_{\epsilon \to 0^+} \int_{\rightarrow} \omega = \int_{1}^{\infty} (-\log x) \biggl( \frac{1}{x-1} - \frac{1}{x+1} \biggr) \, \mathrm{d}x = -\frac{\pi^2}{4}. \end{align*}
Conclusion. Combining all these observations altogether, we have
$$ \bbox[10px,#ffd]{ I = 2 (I_{\rightarrow} + I_{\uparrow} + I_{\swarrow}) = \pi \log 2 } $$
as desired.
${}^{1)}$ We include the sketch of computation for $\mathrm{d}\omega$, which is essentially @G Cab's computation in backward direction. Substituting $u=\theta-\phi$ and $v=\theta+\phi$,
$$ \omega = \left( \log \left( \frac{\sin v - \sin u}{\sin v + \sin u} \right) \right)\biggl( \frac{u \, \mathrm{d}u}{\sin^2 u} - \frac{v \, \mathrm{d}v}{\sin^2 v} \biggr) . $$
Taking $\mathrm{d}$,
\begin{align*} \mathrm{d}\omega &= \mathrm{d}\left( \frac{u}{\sin^2 u} \log \left( \frac{\sin v - \sin u}{\sin u + \sin v} \right) \right) \wedge \mathrm{d}u \\ &\qquad - \mathrm{d}\left( \frac{v}{\sin^2 v} \log \left( \frac{\sin v - \sin u}{\sin u + \sin v} \right) \right) \wedge \mathrm{d}v \\ &= \frac{u}{\sin^2 u} \left( \frac{\partial}{\partial v} \log \left( \frac{\sin v - \sin u}{\sin u + \sin v} \right) \right) \mathrm{d}v \wedge \mathrm{d}u\\ &\qquad - \frac{v}{\sin^2 v} \left( \frac{\partial}{\partial u} \log \left( \frac{\sin v - \sin u}{\sin u + \sin v} \right) \right) \mathrm{d}u \wedge \mathrm{d}v \\ &= 4 \left( \frac{v \cos u \csc v - u \cos v \csc u}{\cos(2u) - \cos (2v)} \right) \, \mathrm{d}u \wedge \mathrm{d}v. \end{align*}
By noting that $\mathrm{d}u \wedge \mathrm{d}v = 2 \, \mathrm{d}\theta \wedge \mathrm{d} \varphi$, the question boils down to showing
$$ 2 \left( \frac{v \cos u \csc v - u \cos v \csc u}{\cos(2u) - \cos (2v)} \right) = \frac{(v+u)\cot(v+u) - (v-u)\cot(v-u)}{\cos(v+u) - \cos(v-u)}. $$
This can be verified by expanding the right-hand side.