Mean distance between 2 points in a ball I have found an answer on this site to the question of determining the mean straight-line distance between 2 randomly chosen points in a disc of radius r. (See Average distance between two points in a circular disk) I'm now trying to find an answer to the same question except involving a ball of radius r rather than a disc. Any guidance on this question would be appreciated.
 A: The average distance can be calculated as a triple integral in polar coordinates:

We ran into it some years back when doing research on proteins (http://www.ncbi.nlm.nih.gov/pubmed/9514112).
A: (a) Two random points on the unit sphere $S^2$:
Assume the first point at the north pole ($\theta=0$) of $S^2$. Then the distance to a point at latitude $\theta$ is $2\sin{\theta\over2}$. Therefore the mean distance between the north pole and the second point is given by
$${1\over 4\pi}\int_0^\pi 2\sin{\theta\over2}\cdot 2\pi\sin\theta\ d\theta={4\over3}\ .$$
$$ $$
(b) Two random points in the unit ball of ${\mathbb R}^3\ $:
Let ${\bf X}$ and ${\bf Y}$ be the two random points. Then $R:=|{\bf X}|$, $\ S:=|{\bf Y}|$, and $\Theta:=\angle({\bf X},{\bf Y})$ are independent random variables with densities
$$f_R(r)=3r^2\quad (0\leq r\leq 1)\ ,\qquad f_S(s)=3s^2\quad(0\leq s\leq 1)\ ,$$
and
$$f_\Theta(\theta)={1\over2}\sin\theta\quad(0\leq\theta\leq\pi)\ .$$
(Concerning $f_R$ and $f_S$ note that the volume included between $r$ and $r+dr$ is proportional to $r^2$. For $f_\Theta$ you may assume ${\bf X}$ pointing due north. The abstract surface area between $\theta$ and $\theta+d\theta$ is then proportional to $\sin\theta$, as in (a).)
It follows that the mean distance $\delta$ between ${\bf X}$ and ${\bf Y}$ is given by
$$\delta=\int_0^1\int_0^1\int_0^\pi \sqrt{r^2+s^2-2rs\cos\theta}\ f_R(r) f_S(s)f_\Theta(\theta) d\theta\ ds\ dr\ .$$
The innermost integral computes to
$$\eqalign{{1\over2}\int_0^\pi \sqrt{r^2+s^2-2rs\cos\theta}\ \sin\theta\ d\theta&={1\over 6rs}\bigl(r^2+s^2-2r s\cos\theta\bigr)^{3/2}\Biggr|_0^\pi \cr 
&={1\over 6rs}\bigl((r+s)^3-|r-s|^3\bigr)\ .
\cr}$$
In the sequel we assume $s\leq r$ and compensate this by a factor of $2$. We are then left with
$$\delta=\int_0^1\int_0^r 3r s(6r^2 s +2s^3)\ ds\ dr={36\over35}\ .$$
It should not be too difficult to set a similar computation up that is valid for a ball in any ${\mathbb R}^n$, $\ n\geq 2$.
A: (a) UNIT SPHERE
Let us recall the Laplace's multipole expansion formula into Legendre polynomials:
$$\frac{1}{|\vec{r}'-\vec{r}|}=\frac{1}{r}\sum_{n=0}^\infty \left(\frac{r'}{r}\right)^n P_n(Y)\qquad ;r\geq r',$$
where I have denoted $Y:=\hat{r}\bullet\hat{r}'$. If we assume that $r=r'=1$, then
$$\frac{1}{|\hat{r}'-\hat{r}|}=\sum_{n=0}^\infty P_n(Y)$$
The average distance between two uniformly distributed points on the unit sphere we denote as $\bar{L}$. Since
$$|\vec{r}'-\vec{r}|=(\vec{r}'-\vec{r})^2/|\vec{r}'-\vec{r}|=(r'^2+r^2-2rr'Y)/|\vec{r}'-\vec{r}|,$$
we have for $r=r'=1$:
$$|\hat{r}'-\hat{r}|=2(1-Y)/|\hat{r}'-\hat{r}|$$
and via Laplace's expansion therefore
$$\bar{L} = \mathbb{E}\left[|\hat{r}'-\hat{r}|\right]=\frac{1}{(4\pi)^2}\oint \oint |\hat{r}'-\hat{r}| d\Omega' d\Omega=\frac{1}{8\pi^2}\sum_{n=0}^\infty\oint \oint \left(1-Y\right) P_n(Y) d\Omega' d\Omega.$$
Due to orthogonality, we know as well that
$$\oint P_n(Y)P_m(Y) d\Omega = \frac{4\pi}{2n+1} \delta_{nm}.$$
Since $1=P_0(Y)$ and $Y = P_1(Y)$, we immediately get
$$ \bar{L} = \frac{1}{8\pi^2}\sum_{n=0}^\infty\oint \oint \left(1-Y\right) P_n(Y) d\Omega' d\Omega = \frac{1}{2\pi} \oint \left(1-\frac{1}{3}\right) d\Omega = \frac{1}{2\pi}\frac{2}{3}4\pi = \frac{4}{3}$$
(b) UNIT BALL
The average distance between uniformly distrubuted two random points in the unit ball is given by$\bar{L} = \mathbb{E}\left[|\vec{r}'-\vec{r}|\right]$. Since the problem is symmetric with respect of labeling of the two points and changing their place, we can define a random variable $L' := |\vec{r}'-\vec{r}|\theta(r-r')$ which vanishes when $r'>r$. Clearly, we miss exactly half of the points, therefore
$$\bar{L} = 2\bar{L}' = 2\mathbb{E}\left[|\vec{r}'-\vec{r}|\theta(r-r')\right]=2\left(\frac{3}{4\pi}\right)^2\int \int |\vec{r}'-\vec{r}|\theta(r-r') dV' dV$$,
which can be expanded by the same trick as before. One gets
$$ \bar{L} = \frac{9}{8\pi^2}\int_0^1 \int_0^r \oint \oint \left(r'^2+r^2-2rr'Y\right) \frac{1}{r} \sum_{n=0}^\infty \left(\frac{r'}{r}\right)^n P_n(Y) r^2r'^2 d\Omega' d\Omega dr' dr$$
Due to orthogonality, we get
$$ \bar{L} = 18\int_0^1 \int_0^r\frac{1}{r}\left(r'^2+r^2-\frac{2}{3}rr'\frac{r'}{r}\right) r^2r'^2 dr'dr = 18\int_0^1 r^6\left(\frac{1}{5}+\frac{1}{3}-\frac{2}{3}\frac{1}{5}\right)dr = \frac{36}{35}$$
(c) BALL DISTRIBUTION FUNCTION
Denote the probablity density of $L$ as a function $g(\lambda)$, for this function one has
$$g(\lambda) = \left(\frac{3}{4\pi}\right)^2\int\int \delta\left(|\vec{r}'-\vec{r}|-\lambda\right)dV'dV.$$
We expand the $\delta(L-\lambda)$ as a sum of Legendre polynomials. Note that $\delta(L-\lambda) \neq 0$ only when $|r'-r|<\lambda<r+r$ since $|\vec{r}'-\vec{r}|\in (|r'-r|,r'+r)$. Since $L = |\vec{r}'-\vec{r}|=\sqrt{r'^2 + r^2 - 2rr'Y}$ and $\delta(f(x))=\delta(x-x_0)/|f'(x_0)|$ with $f(x_0)=0$, we get
$$\delta(L-\lambda) = \frac{\sqrt{r'^2 + r^2 - 2rr' Y}}{rr'}\delta\left(Y-\frac{r'^2+r^2-\lambda^2}{2rr'}\right)=\frac{\lambda}{rr'}\delta\left(Y-\frac{r'^2+r^2-\lambda^2}{2rr'}\right).$$
Without loss of generality, we assume $r'<r$ (symmetry), $r-r'<\lambda<r+r'$ (vanishment), then
$$\delta(L-\lambda) = \sum_{n=0}^\infty A_n P_n(Y)$$
Integrating with $P_m(Y)$ over $d\Omega$ and using the orthogonality relation, we get
$$A_n = \frac{2n+1}{2} \int_{-1}^1 \delta(L - \lambda) P_n(Y) dY = \frac{2n+1}{2} \int_{-1}^1 \frac{\lambda}{rr'}\delta\left(Y-\frac{r'^2+r^2-\lambda^2}{2rr'}\right) P_n(Y) dY,$$
therefore for $|r-r'|<\lambda<r+r'$:
$$\delta(L-\lambda)= \frac{\lambda}{2rr'} \sum_{n=0}^\infty (2n+1)P_n\left(\frac{r'^2+r^2-\lambda^2}{2rr'}\right)P_n(Y).$$
For the probability density therefore, with help of symmetry
$$g(\lambda) = 2\lambda\left(\frac{3}{4\pi}\right)^2 \sum_{n=0}^\infty \underset{0<r-r'<\lambda<r+r'}{\int \int \oint \oint} \frac{2n+1}{2rr'} P_n\left(\frac{r'^2+r^2-\lambda^2}{2rr'}\right)P_n(Y)  r'^2 r^2 d\Omega' d\Omega dr' dr$$
which is due to orthogonality (only $n=0$ survives)
$$g(\lambda) = 9\lambda \int_{\lambda/2}^{1} \int_{\lambda-r}^{r}\!\! r'r\, dr' dr = \frac{9}{2}\lambda^2 \int_{\lambda/2}^{1} \!\! 2r^2\!-\!\lambda r\, dr = \frac{3}{16}\lambda^2 (2-\lambda)^2(4+\lambda)$$
A: This isn't a complete answer, but a start, a rough one at that.
Given two points $p_0$ and $p_1$ find a third point $p_2$ such that $N=(p_0-p_2)\times(p_1-p_2) \neq \varnothing$ Where $\varnothing$ is the null vector. Define $n = \frac{N}{|N|}$ that is normalize the vector $N$, without normalization the distance won't be correct below.
Now we have a plane where $p$ is any point on it defined as:
$$
n\cdot(p-p_2)=0 \\
n\cdot p-n\cdot p_2=0 \\
n\cdot p = n\cdot p_2\\
$$
Recall the general equation of a plane is $ax+by+cz=d$ where $d$ is the distance from the origin to the plane. Expanding out the above equation gives us
$$a=n_x,b=n_y,c=n_z \\
d=n_xp_{2x}+n_yp_{2y}+n_zp_{2z} $$
Radius of the resulting circle in the plane with a sphere of radius $r$ is given as $R = \sqrt{r^2-d^2} $ You may want to negate $n$ if $d$ is negative as that will allow you to have the center of the circle at $(da,db,dc)$
We haven't defined an axes for our plane so its difficult to map from our 3D points onto it. This is where I'm stopping for now. (I might come back and finish this up) I have a feeling this approach is more of a hassle than just computing $p_0-p_1$ and calculating its length. Namely because we have to normalize a vector ($N$) just to get a distance and thus the radius.
Wiki pages I used for this: Plane and Plane-sphere intersection
