Improper integral of hyperbolic function I'm looking for an elementary derivation of the following formula:

$$\int_0^\infty \frac{\sinh(ax)}{e^{bx}-1}dx=\frac{1}{2a}-\frac{\pi}{2b}\cot\frac{a\pi}{b}$$

...only where $|a|\lt b$, of course, to ensure convergence. Does anyone know of any elementary ways of proving this identity?
I'm fairly sure that a nice, elementary antiderivative does not exist, and I haven't been able to get anywhere with any of the "definite integral tricks" that I have learned thus far.

I give a real method that makes use of series and a bunch of special functions and their associated properties so perhaps this is not quite the elementary method you are after.
Let 
$$I = \int^\infty_0 \frac{\sinh (ax)}{e^{bx} - 1} \, dx, \quad |a| < b.$$
Rewriting the hyperbolic sine function term appearing in the numerator in terms of exponentials, since the improper integral converges, we may write
\begin{align*}
I &= \frac{1}{2} \int^\infty_0 \frac{e^{ax}}{e^{bx} - 1} \, dx - \frac{1}{2} \int^\infty_0 \frac{e^{-ax}}{e^{bx} - 1} \, dx\\
&= \frac{1}{2} \int^\infty_0 \frac{e^{-(b - a)x}}{1 - e^{-bx}} - \frac{1}{2} \int^\infty_0 \frac{e^{-(b + a) x}}{1 - e^{-bx}} \, dx.
\end{align*}  
Recognising the term $\dfrac{1}{1 - e^{-bx}}$ as the sum of a convergent geometric series, namely
$$\frac{1}{1 - e^{-bx}} = \sum^\infty_{n = 0} e^{-bnx}$$
the integral appearing above, after interchanging the summation with the integration, becomes
$$I = \frac{1}{2} \sum^\infty_{n = 0} \int^\infty_0 e^{-(b + nb - a)x} \, dx - \frac{1}{2} \sum^\infty_{n = 0} \int^\infty_0 e^{-(b + nb + a) x} \, dx.$$
Note that as $|a| < b$ the exponents appearing in each of the exponential terms found in each of the integrals will be negative ensuring convergence of the improper integrals. After performing the integrations we are left with
\begin{align*}
I &= \frac{1}{2} \sum^\infty_{n = 0} \frac{1}{(b - a) + nb} - \frac{1}{2} \sum^\infty_{n = 0} \frac{1}{(a + b) + nb}\\
&= \frac{1}{2b} \sum^\infty_{n = 0} \frac{1}{\left (\frac{b - a}{b} \right ) + n} - \frac{1}{2b} \sum^\infty_{n = 0} \frac{1}{\left (\frac{a + b}{b} \right ) + n}.
\end{align*}
Each of these sums can be expressed in terms of the Hurwitz zeta function $\zeta(s,q)$ which is defined as
$$\zeta(s,q) = \sum^\infty_{n = 0} \frac{1}{(q + n)^s}.$$
Thus
$$I = \frac{1}{2b} \zeta \left (1,1 - \frac{a}{b} \right ) - \frac{1}{2b} \zeta \left (1,1 + \frac{a}{b} \right ).$$
Now the Hurwitz zeta function is related to the polygamma function function of order $m$, which is denoted by $\psi^{(m)} (z)$, by (see here)
$$\psi^{(m)} (z) = (-1)^{m + 1} m! \zeta (m + 1, z).$$
Setting $m = 0$ gives
$$\zeta (1,z) = - \psi (z),$$
where $\psi (z) = \psi^{(0)} (z)$ is the digamma function. So in terms of digamma functions our integral can be written as
$$I = \frac{1}{2b} \left [\psi \left (1 + \frac{a}{b} \right ) - \psi \left (1 - \frac{a}{b} \right ) \right ].$$
Since (see here)
$$\psi (z + 1) = \psi (z) + \frac{1}{z},$$
we can write
$$\psi \left (1 + \frac{a}{b} \right ) = \psi \left (\frac{a}{b} \right ) + \frac{b}{a},$$
and from the reflection formula for the digamma function, namely
$$\psi (1 - z) - \psi (z) = \pi \cot (\pi z),$$
we have
$$\psi \left (1 - \frac{a}{b} \right ) = \psi \left (\frac{a}{b} \right ) + \pi \cot \left (\frac{a \pi}{b} \right ),$$
yielding
$$\int^\infty_0 \frac{\sinh (ax)}{e^{bx} - 1} \, dx = \frac{1}{2a} - \frac{\pi}{2b} \cot \left (\frac{a \pi}{b} \right ),$$
as claimed.
