Suppose that $w=\sqrt{R^2 + s^2 -2Rs\cos\theta}$ with $R\ge s>0$, that $p$ is a positive integer and that $j$ is an integer with $0\le j\le p$. Let $I$ and $K$ denote the modified Bessel functions of the first and second kind respectively. I want to prove the following identity (actually I would be happy with the case $j=0$). \begin{align*} \sqrt{\frac{2}{\pi}}&\int_0^\pi \frac{K_{j-1/2}(w)}{w^{j-1/2}}\sin^{2p-1}(\theta)\, d\theta \\ &= (-1)^{p-j}(p-1)! 2^p %\\&\qquad \sum_{m=0}^{p-j}(-1)^m \binom{p-j}{m} R^{2m} \frac{K_{p+m-1/2} (R)}{R^{p+m-1/2}}\cdot \frac{I_{m+j-1/2} (s)} { s^{m+j-1/2}} \end{align*} This is very reminiscent of expressions you see in the chapters on addition formulas and definite integrals in Waton's book on Bessel functions, but I have no idea about how to prove it; I have, however, verified it numerically with Sage for small values of $p$ and $j$.
The case that $j=0$ reduces to \begin{align*}\int_0^\pi &e^{-w}\sin^{2p-1}(\theta)\, d\theta \\&= (-1)^{p}(p-1)! 2^p \sum_{m=0}^{p}(-1)^m \binom{p}{m} R^{2m} \frac{K_{p+m-1/2} (R)}{R^{p+m-1/2}}\cdot \frac{I_{m-1/2} (s)} { s^{m-1/2}}, \end{align*} and this is equivalent to the identity that I was asking about in Proof of an identity involving $\exp(−|x−s|)dx$ over an even sphere. I'm hoping that my generalization to the above expression, which looks like a typical Bessel function identity, will be more familiar to people and be open to more standard methods of proof.