Integrate $\frac{R}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} \frac{1-\cos(mx+ny)}{2-(\cos x + \cos y)} dx dy $ As in the title: let $m,n\in\mathbb{Z}$. Integrate:
$$ \frac{R}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} \frac{1-\cos(mx+ny)}{2-(\cos x + \cos y)} dx dy $$ 
For me, it's quite a difficult problem. Any hints? 
$R$ is a constant. 
 A: Approach 
I will show that when $R=1$ the resulting discrete scalar distribution:
$$ A(m,n) = \frac{1}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} \frac{1-\cos(mx+ny)}{2-(\cos x + \cos y)} dx dy $$
is harmonic everywhere except when $(m,n)=(0,0)$. (In fact, $A (m,n)$ is the inverse  operator of the Discrete Laplacian).
This allows us to find all integrals with recursion, once the simple forms are calculated. 
The discrete Laplacian 
Let: 
$$ g(m,n) = \frac{1-\cos(mx+ny)}{2-(\cos x + \cos y)} $$
Then the Discrete Laplacian of $g(m,n)$  in the $(m,n)$ domain is:  
$$ \begin{align} \mathbf{D}^2_{mn} \otimes g(m,n) &=\begin{bmatrix}0 & \frac 1 2 & 0\\ \frac 1 2 & -2 & \frac 1 2\\0 & \frac 1 2 & 0\end{bmatrix} \otimes g(m,n) \\ & = \cos(mx+ny) \end{align}$$
This is not to hard to proof: Wolfram alpha

Harmonicity 
Now define:
$$ \begin{align} I(m,n) &\equiv \mathbf{D}^2_{mn} \otimes A (m,n) \\&=\mathbf{D}^2_{mn} \otimes \frac{1}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} g(m,n)dx dy \\ & =  \frac{1}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} \cos(mx+ny) dx dy \\ &= \mathcal{F}(1)\\ &= {\begin{cases}
1 & \text {if $(m,n)=(0,0)$} \\
0 & \text{otherwise} \\
\end{cases}} \end{align} $$
which is the identity operator for convolution. This means that for each pair of integers $(m,n) \neq (0,0)$ the integral is the average of the four integrals of its direct neighbors $\{(m-1,n), (m+1,n), (m,n-1), (m,n+1)\}$. 

Startvalues
It is easy to see that for $(m,n) = (0,0)$ the integral is zero. For the four direct neigbors $m^2 + n^2 = 1$, the integral must be $\frac 12$, because the Laplacian at $(0,0)$ is one, and because of symmetry.   
Together with the integral derived for $m=m$ in the
answer of Sangchul we have enough initial values to recursively find the other integrals:
$$ \frac{1}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} g(m,n) dx dy = {\begin{cases}
0 & \text {if $(m,n)=(0,0)$} \\
0.5 & \text {if $ m^2+n^2 = 1$} \\
\sum_{k=1}^m  \frac {2} {\pi(2 k - 1)} & \text {if $ m=n$} \\ \end{cases}} $$
Recursive recipe:
If $n \neq m$, assume that $n>m$ (symmetry) then: 
$$A(m,n) =  4*A(m,n-1) - (A(m-1,n-1) + A(m+1,n-1) + A(m,n-2)) $$ 
$\blacksquare$
A: Just an addendum to Sangchul Lee's answer.
If we exploit the parity of the cosine function we get that the wanted integral equals:
$$ I(m,n) = \frac{R}{\pi^2}\int_{0}^{\pi}\int_{0}^{\pi}\frac{1-\cos(nx)\cos(my)}{2-(\cos x+\cos y)}\,dy\,dy $$
and by writing $\frac{1}{2-(\cos x+\cos y)}$ as $\int_{0}^{+\infty}\exp\left[\left(\cos x+\cos y-2\right)t\right]\,dt$ we get:
$$ I(m,n) = R\int_{0}^{\infty}e^{-2t}\left(I_0(t)^2-I_n(t)\,I_m(t)\right)\,dt $$
where $I_k$ is a modified Bessel function of the first kind. Since the inverse Laplace transform of $e^{-2t}$ is just $\delta(s-2)$, the problem boils down to computing the Laplace transform of $I_0(t)^2$ and $I_n(t)\,I_m(t)$. The Laplace transform of $I_0(t)^2$ is related with the complete elliptic integral of the first kind, and has a singularity at $s=2$, where it behaves like $-\frac{\log|2-s|}{2\pi}$.
In a explicit way, from:
$$ I_0(z) = \sum_{n\geq 0}\frac{z^{2n}}{4^n n!^2} $$
it follows that, for any $\tau>2$:
$$ I_0(z)^2 = \sum_{n\geq 0}\frac{\binom{2n}{n}}{4^n n!^2}\,z^{2n},\qquad \int_{0}^{+\infty}I_0(z)^2 e^{-\tau z}\,dz = \sum_{n\geq 0}\frac{\binom{2n}{n}^2}{4^n \tau^{2n+1}}=\frac{2}{\pi\tau}\,K\left(\frac{2}{\tau}\right) $$
At least in principle, we may perform the same computation for $I_m(t)\,I_n(t)$, then compute the Laplace transform of $I_0(t)^2-I_m(t)\,I_n(t)$ a $\tau=2$ through De l'Hopital theorem.
