An integral from Newton's theorem One way of proving Newton's theorem on the potential of a spherically symmetric measure (see below) is to evaluate the integral
$$ \int_0^{\pi} \frac{\sin^{n-2}{\theta}}{(1+r^2 - 2r \cos{\theta} )^{(n-2)/2}} \, d\theta = \frac{\sqrt{\pi} \Gamma((n-1)/2)}{\Gamma(n/2)} (\max{(1,r)})^{2-n}. \tag{1}$$
In what ways can one prove this for general $n \neq 2$? (Integer $n$ is obviously the interesting case.) I'm looking in particular for some more elementary ways of doing it (I've added a couple of rather high-powered ways below), but otherwise, anything goes: complex analysis, Fourier transforms, whatever.
(The case $n=2$ is much easier, and has been answered before on this site.)
It is clear that for $r>1$, we can take out a factor of $r^{2-n}$, reducing to the $r<1$ case ($r=1$ is easy: just apply double-angle formulae and the definition of the beta function).

(The following is the context from which the above arises.)
Let $\mu(x)$ be a measure on $\mathbb{R}^n$. The Coulomb potential of $\mu$, $D(\mu)(x)$, is defined as
$$ D(\mu)(x) = \int_{\mathbb{R}^n} G_n(x-y) \, d\mu(y), $$
where $G_n$ is the fundamental solution of Laplace's equation in $\mathbb{R}^n$,
$$ G_n(x) = \begin{cases} \frac{1}{(n-2)S_{n-1}} \lvert x \rvert^{2-n} & n \neq 2 \\
\frac{1}{2\pi} \log{\lvert x \rvert} & n=2
\end{cases}, $$
$S_{n-1}$ the surface area of the $(n-1)$-dimensional sphere.
(The point being of course that $-\Delta D(\mu) = \mu $.)
Newton's theorem states that for a spherically symmetric measure $\mu$, the Coulomb potential satisfies
$$ D(\mu)(x) = G_n(x) \int_{|y|\leqslant |x|} d\mu(y) + \int_{|y|>|x|} G_n(y) \, d\mu(y) $$
(see, for example, Lieb and Seiringer, The Stability of Matter in Quantum Mechanics, p.91f.) (We assume that all the appropriate integrals converge.)
Now, one way of proving this (indeed, probably the simplest way) is to turn the integral into spherical coordinates and do the integrals over spheres. It is easy to see that $\lvert x-y \rvert = \sqrt{\lvert x \rvert^2 + \lvert y \rvert^2 - 2\lvert x \rvert \lvert y \rvert \cos{\theta} }$, $\theta$ the angle between $x$ and $y$. Choosing the right spherical coordinates eventually gives us the integral
$$ S_{n-2}\int_0^{\pi} \frac{\sin^{n-2}{\theta}}{(\lvert x \rvert^2 + \lvert y \rvert^2 - 2\lvert x \rvert \lvert y \rvert \cos{\theta} )^{(n-2)/2}} \, d\theta, $$
for $n \neq 2$, which we need to show is equal to
$$ S_{n-1} \max{(\lvert x\rvert ,\lvert y\rvert )})^{2-n}, $$
and reducing further leads to (1).
 A: One can employ Sturm-Liouville theory and Gegenbauer polynomials: they have generating function
$$ \frac{1}{(1-2xt+t^2)^{\alpha}} = \sum_{n=0}^{\infty} t^n C_{n}^{(\alpha)}(x), $$
and form an orthogonal system of polynomials with weight $(1-x^2)^{\alpha-1/2}$. Hence if $|x|<1$,
$$ \int_{-1}^1 \frac{(1-x^2)^{\alpha-1/2}}{(1-2xt+t^2)^{\alpha}} \, dx = \sum_{n=0}^{\infty} \int_{-1}^1 C_{0}^{(\alpha)}(x) C_{n}^{(\alpha)}(x) (1-x^2)^{\alpha-1/2} \, dx = \int_{-1}^1 (1-x^2)^{\alpha-1/2} \, dx = \frac{\sqrt{\pi}\Gamma(\alpha+1/2)}{\Gamma(\alpha+1)}, $$
which is what we want, taking $\alpha=n/2-1$: this is the same as our integral after setting $x=\cos{\theta}$.

Alternatively, the transformation $u=\frac{1}{2}(1+\cos{\theta})$ puts the integral in the form
$$ 2^{n-2}\int_0^1 (u(1-u))^{(n-3)/2} \left((1+r)^2-4ru \right)^{-(n-2)/2} \, du, $$
which is almost a hypergeometric function,
$$ 2^{n-2} (1+r)^{2-n} \frac{\Gamma((n-1)/2)^2}{\Gamma(n-1)} F\left( {n/2-1,n/2-1/2 \atop n-1} ; \frac{4r}{(1+r)^2} \right) $$
We have (DLMF 15.8.E15) the quadratic transformation
$$ F\left({a,b\atop a-b+1};z\right)=(1+z)^{-a} F \left({\frac{1}{2}a,\frac{1}{2}a+\frac{1}{2}\atop a-b+1};\frac{4z}{%
(1+z)^{2}}\right),$$ 
which we can use with $a=n-2,b=0$. Then the left-hand side of the identity is just $1$, and using some Gamma-function identities gives the correct answer: $\Gamma((n-1)/2)/\Gamma(n-1) = 2^{2-n}\sqrt{\pi}/\Gamma(n/2)$.
A: I don't have time to write a full explanation here, but (unless I misunderstand your question) a simple argument is outlined in Section $9.7$ of the book of Lieb and Loss (Analysis, 2nd edition, AMS, Grad. Stud. in Math., vol. 14), based on the observation that for fixed $y\in\mathbb{R}^N$ the function $$x\mapsto |x-y|^{-(n-2)}$$ is harmonic on $\{x:|x|<|y|\}$ (and the ensuing mean value property).
