Roll a sphere over two other spheres which it contacts Problem statement:
There are three spheres. The one which will roll is $\textbf{X}=(x_1,x_2,x_3)$ with radius $R_X$. The other two spheres are $\textbf{A}=(a_1,a_2,a_3)$ with $R_A$ and $\textbf{B}=(b_1,b_2,b_3)$ with $R_B$. They are both below $X$, meaning $x_3>a_3$ and $x_3>b_3$. They both contact $X$, meaning $|\textbf{X}-\textbf{A}| = R_X+R_A$ and similarly for $B$. I need to find the axis of rotation $\textbf{V}$ and point $\textbf{v}$ about which $\textbf{X}$ can rotate to change its position while preserving its contacts with both $\textbf{A}$ and $\textbf{B}$. With this information actually performing the rotation for my task in a granular physics program is simple. 
Attempted solution:
I thought the axis of rotation would be the vector connecting $A$ and $B$. 
$$\textbf{V}= \textbf{B}-\textbf{A},$$
suitably normalized. I thought the point about which the rotation should occur would come from analyzing the triangle formed by the locations of the three spheres, and breaking it into two right triangles with height $h$ and bases $d$ and $D$, then solving the system
$$(R_A+R_X)^2 = d^2 + h^2,$$
$$(R_B + R_X)^2 = D^2 + h^2,$$
$$d+D = |\textbf{B}-\textbf{A}|,$$
which is three equations in three unknowns $h$,$d$, and $D$. I thought $h$ would be the length of the lever arm about which to rotate, and I thought $d$ or $D$ would map  $\textbf{A}$ to the pivot point as $\textbf{v} = \textbf{V}d$. When I use this for the rotation I get nonsense results. What am I missing? Any help is appreciated!  
Here is a figure to indicate a possible configuration before rolling. It is not necessary however that the two lower spheres $\textbf{A}$ and $\textbf{B}$ are touching. 
It works. Thanks!

 A: By translating coordinates we may as well take ${\bf B} = {\bf 0}$, so that $\ell = \langle {\bf A} \rangle$. Then, there is a (unique) orthogonal decomposition $${\bf X} = {\bf X}^{\parallel} + {\bf X}^{\perp} ,$$
where ${\bf X}^{\parallel} \parallel \ell$ and ${\bf X}^{\perp} \perp \ell$. (Explicitly, ${\bf X}^{\parallel} = \textrm{proj}_{\bf A} {\bf X}$.)
Now, as I pointed out in my comment, any rotation fixing the spheres $A$ and $B$ fixes ${\bf A}, {\bf B}$, so it also fixes the unique line $\ell$ through them and hence fixes ${\bf X}^{\parallel}$. On the other hand, rotations about $\ell$ also fix the plane $\Pi$ orthogonal to $\ell$ and $\bf X$, which by construction intersects $\ell$ at ${\bf X}^{\parallel}\textrm{.}$ The "circular track" along which ${\bf X}$ moves under rotation mentioned in the comments is the circle in $\Pi$ centered at $||{\bf X}^{\parallel}||$ of radius $||{\bf X}^{\perp}||$ (this circle degenerates to a single point in the case that ${\bf X} \in \ell$, that is, when the three centers are collinear.
A: The axis of rotation is, as you’ve already determined, the line through the center of the two fixed spheres, which can be represented parametrically as $(1-\lambda)\mathbf A+\lambda\mathbf B$. This line has direction vector $\mathbf B-\mathbf A$, which is normalized to $\mathbf w=\|\mathbf B-\mathbf A\|$. When a point that’s not on this axis is rotated about the axis, it moves along a circular arc that lies in a plane perpendicular to the rotation axis and is centered on that axis.  
To find this rotation center $\mathbf C$ for the center of the moving sphere $\mathbf X$, observe from the triangle of the sphere centers that $\mathbf C$ is the orthogonal projection of $\mathbf X$ onto $\overline{\mathbf A\mathbf B}$. Therefore, $$\mathbf C=\mathbf A+\mathbf w\cdot(\mathbf X-\mathbf A)\mathbf w.$$ Here, $\mathbf w\cdot(\mathbf X-\mathbf A)$ is the distance $d$ that you were solving for. This point could also be found by directly computing the intersection of the rotation axis and its perpendicular plane through $\mathbf X$. One way to do this is to solve the equation $$[(1-\lambda)\mathbf A+\lambda\mathbf B]\cdot\mathbf w=\mathbf X\cdot\mathbf w$$ for $\lambda$, which leads directly to the above result. (This equation comes from the point-normal form of equation for a plane.)  
It seems useful to actually find the circle of motion. Let $\mathbf u=\|\mathbf X-\mathbf C\|$ and $\mathbf v=\mathbf w\times\mathbf u$. These two unit vectors are both parallel to the plane of motion and are orthogonal to each other. A parametric description of $\mathbf X$’s motion is then $\mathbf C+R_X\mathbf u\cos\theta+R_X\mathbf v\sin\theta$, a clockwise rotation through $\theta$ when viewed from $\mathbf A$. (I didn’t really need to normalize $\mathbf u$ for this, but I’m anticipating the next paragraph.)  
With the center point $\mathbf C$ and the three unit vectors in hand, it’s fairly easy to construct a homogeneous rotation matrix that implements this rotation. In the $uvw$ coordinate system with origin at $\mathbf C$, it’s just a rotation about the $w$-axis, which can be represented in block form as $\small\begin{bmatrix}R_\theta&\mathbf0_2\\\mathbf0_2&I_2\end{bmatrix}$, where $R_\theta$ is the familiar two-dimensional rotation matrix and $\mathbf0_2$ is a $2\times2$ matrix of zeros. The corresponding matrix in the $xyz$ coordinate system is obtained via a change of basis: $$R=\begin{bmatrix}\mathbf u&\mathbf v&\mathbf w&\mathbf C \\ 0&0&0&1\end{bmatrix}\begin{bmatrix}R_\theta&\mathbf0_2\\\mathbf0_2&I_2\end{bmatrix}\begin{bmatrix}\mathbf u&\mathbf v&\mathbf w&\mathbf C \\ 0&0&0&1\end{bmatrix}^{-1}.$$ The change-of-basis matrix is a composition of a rotation and translation, so its inverse can be computed without going through a generic matrix inversion: $$R=\begin{bmatrix}\mathbf u&\mathbf v&\mathbf w&\mathbf C \\ 0&0&0&1\end{bmatrix}\begin{bmatrix}R_\theta&\mathbf0_2\\\mathbf0_2&I_2\end{bmatrix}\begin{bmatrix}\mathbf u^T&0\\\mathbf v^T&0\\\mathbf w^T&0\\\mathbf0^T&1\end{bmatrix}\begin{bmatrix}I_3&-\mathbf C\\\mathbf0&1\end{bmatrix}.$$ N.B.: in the first matrix, the unit vectors are columns, while in the third they are rows; $\mathbf C$ and $-\mathbf C$ are both entered as columns. Note, too, that any point on the rotation axis can be used instead of $\mathbf C$. You will end up with the same matrix. (If you like, expand the first matrix into its rotation and translation components to restore the symmetry of the product.)
