Finding circle of a sphere through two points We have two points $P_1, P_2$ on a sphere $S$ of radius $R$. Suppose for $r \ll R$, the distance between $P_1$ and $P_2$ is less than $2r$. 
Then, $P_1$ and $P_2$ both lie on exactly two  radius-$r$ circles of the sphere. This is clear: of the set of planes that contain $P_1$ and $P_2$, there are two of those planes that intersect $S$ to form a circle with radius $r$. 
Given $P_1$, $P_2$, $R$ and $r$, how can I calculate the centers of these two circles?
I would prefer to do as little messy spherical geometry as possible :).

I'm assuming that you want a circle in space, with its center inside $S$ and radius $r$ measured along a straight line in space. If you instead want a circle on the sphere, i.e. with center a point on the sphere, and radius measured along a geodesic arc, this answer doesn't exactly apply. But for $r\ll R$ it might still be useful. Otherwise, projecting the center in space back onto the sphere is easy, so all you'd have to do is translate the geodesic radius to a straight line radius.
The distance between the center of the sphere and the center of one of these circles has to be
$$ d = \sqrt{R^2 - r^2} $$
This can be seen from the following cross section through the sphere:

So you are looking for a point $C$ such that $\lVert P_1-C\rVert=\lVert P_2-C\rVert=r$ and $\lVert C\rVert = d$. If we use coordinates
$$P_1=(x_1,y_1,z_1)\qquad P_2=(x_2,y_2,z_2)\qquad C=(x,y,z)$$
then you can rewrite this to
$$(x_1-x)^2 + (y_1-y)^2 + (z_1-z)^2 = (x_2-x)^2 + (y_2-y)^2 + (z_2-z)^2 = r^2\\
x^2 + y^2 + z^2 = R^2 - r^2$$
At that point you can fire up your computer algebra system to find the solutions here. If you want to continue manually, rewrite the equatoons like this:
\begin{align*}
x_i^2 - 2x_ix + x^2 + y_i^2 - 2y_iy + y^2 + z_i^2 - 2z_iz + z^2 &= r^2 \\
r^2 + 2x_ix + 2y_iy + 2z_iz - x_i^2 - y_i^2 - z_i^2 &= x^2 + y^2 + z^2 \\
r^2 + 2x_ix + 2y_iy + 2z_iz - R^2 &= R^2 - r^2 \\
x_ix + y_iy + z_iz &= R^2 - r^2
\end{align*}
So at that point you have two linear equations:
$$\begin{pmatrix}
x_1 & y_1 & z_1 \\
x_2 & y_2 & z_2
\end{pmatrix}\cdot C=\begin{pmatrix}
R^2 - r^2 \\
R^2 - r^2
\end{pmatrix}$$
Since you have two equations for three variables, you get a one-dimensional space of solutions: $C\in\{v_1+\lambda v_2\mid\lambda\in\mathbb R\}$. Now revisit the last equation:
$$x^2+y^2+z^2 = \lVert C\rVert^2 = \langle v_1+\lambda v_2,v_1+\lambda v_2\rangle = \langle v_1,v_1\rangle + 2\lambda\langle v_1,v_2\rangle + \lambda^2\langle v_2,v_2\rangle\\= \lVert v_1\rVert^2 + 2\lambda\langle v_1,v_2\rangle + \lambda^2\lVert v_2\rVert^2 = d^2 = R^2 - r^2$$
So that's a quadratic equation in $\lambda$, and its two solutions describe your two circle centers.
