In this paper about point vortices in a doubly periodic domain, the authors encounter the following double sum in equation (4), namely $$ S_{M,N}(x,y) = \sum_{m=-M}^M \sum_{n=-N}^N \frac{x-2\pi n}{(x-2\pi n)^2+(y-2\pi m)^2} \tag{4} $$ and $S(x,y) = \lim_{N,M\to \infty} S_{M,N}(x,y)$.
In the words of the authors of the paper, "Using the method of Laplace transforms on the sum over $n$, treating resulting integrals involving $M$ and $N$ asymptotically, and then taking the limit $M, \, N \to \infty$ only for those terms that are independent of the order of the limit results in
$$ S_{M,N} = \sum_{m=-\infty}^\infty \frac{1}{2} \frac{\sin(x)}{\cosh(y-2\pi m)-\cos(x)}+ \frac{x}{2\pi} - \frac{x}{\pi^2}\arctan(N/M). " \tag{5} $$
I am struggling to reproduce this result and I would appreciate any hint that could point me in a fruitful direction. I am not terribly interested in justifying interchanging integrals and infinite sums, I just want to see how to calculation works roughly, if necessary I can fill in the gaps by myself. Here is what I have tried so far.
I know the method of Laplace transforms for summing infinite series: $$ \sum_{n=-N}^N \frac{x-2\pi n}{(x-2\pi n)^2+(y-2\pi m)^2} = \\ -\frac{x}{x^2+(y-2\pi m)^2} +\sum\limits_{n=0}^N \underbrace{\left[ \frac{x-2\pi n}{(x-2\pi n)^2+(y-2\pi m)^2} + \frac{x+2\pi n}{(x+2\pi n)^2+(y-2\pi m)^2} \right]}_{\equiv g(n)} $$
The inverse Laplace transform of $g(n)$ may be calculated to be: $$ G(k)=(\mathcal{L}^{-1}_n[g(n)])(k) =- \frac{1}{\pi} \sinh\left(\frac{xk}{2\pi} \right) \cos\left[\left(\frac{y-2\pi m}{2\pi}\right)k\right]$$
Hence, using that $ g(n) = \int_0^\infty G(n) e^{-kn}\mathrm{d}k $, and the partial sum formula of the geometric series, we get:
$$S_{M,N} = \sum_{m=-M}^M \left[ -\frac{x}{x^2+(y-2\pi m)^2}-\frac{1}{\pi} \int_0^\infty \frac{(1-e^{-(N+1)k})}{1-e^{-k}} \sinh\left(\frac{xk}{2\pi} \right) \cos\left[\left(\frac{y-2\pi m}{2\pi}\right)k\right]\mathrm{d}k \right] $$
I can calculate the infinite sum over the first term in the limit $M\to\infty$, finding $$ \sum_{m=-\infty}^\infty -\frac{x}{x^2+(y-2\pi m)^2} = \frac{\sinh(x)}{2(\cos(y)-\cosh(x))}. $$ Thus, I am at
$$ S_{M,N} \sim \frac{\sinh(x)}{\cos(y)-\cosh(x)} + \sum_{m=-M}^M \int_0^\infty \frac{(1-e^{-(N+1)k})}{1-e^{-k}} \sinh\left(\frac{xk}{2\pi} \right) \cos\left[\left(\frac{y-2\pi m}{2\pi}\right)k\right]\mathrm{d}k$$
Now, it seems tempting to pull the sum over $m$ into the integral and use
$$ \sum_{m=-M}^M \cos\left[\left(\frac{y}{2\pi}-m\right)k\right] = \csc\left( k/2\right) \sin\left( k(M+1/2) \right) \cos\left( \frac{ky}{2\pi}\right)$$
Now I could evaluate the remaining integral as a stationary-phase type integral since the term $ \sin\left( k(M+1/2) \right)$ is rapidly oscillating as $M\to \infty$. However, I don't think this is how (5) was found since it still contains a sum over $m$. Does anyone have an idea how to proceed here?
P.S.: The authors reference a book "A. D. Wheelon, Tables of Summable Series and Integrals Involving Bessel Functions", but that only describes the general Laplace Transform method, not the specific problem at hand.