Supose we have the 3-ball of radii 2, $\mathbb{B}^3_{r=2}=\{ \mathbf{x}\in\mathbb{R}^3: |\mathbf{x}| \leq 2 \}$, and the 2-sphere $\mathbb{S}^2=\{ \mathbf{x}\in\mathbb{R}^3: |\mathbf{x}| = 1 \}$. We know that $\{ \mathbf{w}=\mathbf{v} - \mathbf{u}: \mathbf{v}\in \mathbb{S}^2,\mathbf{u}\in \mathbb{S}^2 \} = \mathbb{B}^3_{r=2}$ thus every point on $\mathbb{B}^3_{r=2}$ can be decomposed into two point on $\mathbb{S}^2$.
Let us say we introduce a parameterization of $\mathbb{S}^2$, $\mathbf{k}(\theta,\phi) = \begin{pmatrix} \cos(\theta)\sin(\phi) \\ \sin(\theta)\sin(\phi) \\ \cos(\phi) \end{pmatrix}$ so that $\{\mathbf{k}(\theta,\phi) \in \mathbb{R}^3 \;\forall\; (\theta,\phi) \in [0,2\pi] \times [0,\pi]\} = \mathbb{S}^2$.
Given a point $\mathbf{b}_0 \in \mathbb{B}^3_{r=2}$, how do we find all combinations of $(\theta_1,\phi_1)$ and $(\theta_2,\phi_2)$ so that $\mathbf{k}(\theta_1,\phi_1) - \mathbf{k}(\theta_2,\phi_2) = \mathbf{b}_0$?
My first try was simply solving the equations as (3 equations, 4 variables, set 1 variable as free): $$ \phi_2 = t $$ $$ \cos(\phi_1) - \cos(\phi_2) = b_{0z} \Leftrightarrow \phi_1 = \cos^{-1}(b_{0z} + \cos(t)) $$ $$ \cos(\theta_1)\sin(\phi_1(t)) - \cos(\theta_2)\sin(t) = b_{0x} \Leftrightarrow \theta_1 = \cos^{-1} \left ( \frac{b_{0x} + \cos(\theta_2)\sin(t)}{\sin(\phi_1(t))} \right ) $$ $$ \sin(\theta_1)\sin(\phi_1(t)) - \sin(\theta_2)\sin(t) = b_{0y} \Leftrightarrow \\ \Leftrightarrow \sin(\cos^{-1} \left ( \frac{b_{0x} + \cos(\theta_2)\sin(t)}{\sin(\phi_1(t))} \right ))\sin(\phi_1(t)) - \sin(\theta_2)\sin(t) = b_{0y} \Leftrightarrow \\ \Leftrightarrow \sqrt{ 1 - \left ( \frac{b_{0x} + \cos(\theta_2)\sin(t)}{\sin(\phi_1(t))} \right )^2}\sin(\phi_1(t)) - \sin(\theta_2)\sin(t) = b_{0y} \Leftrightarrow \\ \Leftrightarrow \sin^2(\phi_1(t)) - \left ( \frac{b_{0x} + \cos(\theta_2)\sin(t)}{\sin(\phi_1(t))} \right )^2\sin^2(\phi_1(t)) = (b_{0y}+ \sin(\theta_2)\sin(t))^2 \Leftrightarrow \\ \Leftrightarrow ... \Leftrightarrow \\ \Leftrightarrow \theta_2 = \sin^{-1} \left (\frac{ \sin^2(\phi_1(t)) - b_{0y}^2 -b_{0x}^2 -\sin^2(t)}{2\sqrt{b_{0y}^2 + b_{0x}^2}\sin(t) } \right ) - atan2(b_{0y},b_{0x})$$
And then i have a closed form expression for $\theta_1(t),\phi_1(t),\theta_2(t),\phi_2(t)$ but when i try this formula numerically it does not work as intended. Can anyone identify either my mistake or give a path (reference to a solution or the solution itself) to the answer?
Surly this question has been solved somewhere before but i cannot find it in any of the books on my shelf nor trough google :P