Constructing tangent spheres with centers located on vertices of an irregular tetrahedron

I'm looking to calculate the positions of a cluster of pairwise tangent spheres. As such one can understand that the center of each sphere $S_i$ in the cluster will exist at a vertex of a corresponding irregular tetrahedron with coordinates $(x_i, y_i, z_i)$ in $\mathbb{R}^3$, and will have some radius $R_i$ such that the set of spheres is: $$\begin{cases} (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = R_0^2\\ (x - x_1)^2 + (y - y_1)^2 + (z - z_1)^2 = R_1^2\\ (x - x_2)^2 + (y - y_2)^2 + (z - z_2)^2 = R_2^2\\ (x - x_3)^2 + (y - y_3)^2 + (z - z_3)^2 = R_3^2\\ \end{cases}$$

All sphere radii are known, along with the coordinates of $S_1, S_2, S_3$: the only unknowns are $(x_0,y_0,z_0)$. Since the spheres must be tangent, we also know all side lengths of the tetrahedron: $d(S_i,S_j) = R_i+R_j$.

Solving this set of equations directly has got a little out of hand for me, so I have also attempted to identify the coordinates of $S_3$ by solving the tetrahedron EFGH in a simple basis following this answer.

Once identified, I then attempt to transform them to the original basis ABCD first by

1. Rotating EF onto AB

2. Rotating G to C by rotating face normals of the triangle obtained in 1

3. Translating the result by A

This tends to have a number of edge cases though, for example EFG=ACB has a mirror symmetry that must be checked for.

I'd much prefer a more direct approach (if it exists), similar to this solution without the simplification that the base triangle is equilateral.

Can someone shed some light on this problem for me?

Consider the system \begin{align} (x-x_1)^2+(y-y_1)^2+(z-z_1)^2=d_1^2\tag{1}\\ (x-x_2)^2+(y-y_2)^2+(z-z_2)^2=d_2^2\tag{2}\\ (x-x_3)^2+(y-y_3)^2+(z-z_3)^2=d_3^2\tag{3}\\ \end{align} as a solution to the problem, where $(x,y,z)$ are the unknowns to solve.

Subtract (1) from (2)&(3) to obtain the equations for two planes, and couple this set with (1) to yield the quadratic system: \begin{align} 2(x_1-x_2)x+2(y_1-y_2)y+2(z_1-z_2)z&=d_2^2-d_1^2+x_1^2-x_2^2+y_1^2-y_2^2+z_1^2-z_2^2\\ 2(x_1-x_3)x+2(y_1-y_3)y+2(z_1-z_3)z&=d_3^2-d_1^2+x_1^2-x_3^2+y_1^2-y_3^2+z_1^2-z_3^2\\ x^2+y^2+z^2-2 x_1 x-2 y_1 y-2 z_1 z &= d_1^2-x_1^2-y_1^2-z_1^2\\ \end{align}

Which now must be solved, although the solution seems to blow up in $x$ when I attempt this. We can simplify by eliminating the $S_1$ components such that \begin{align} x^2+y^2+z^2-2x_2 x-2y_2 y-2z_2 z&=d_2^2-x_2^2-y_2^2-z_2^2\\ x^2+y^2+z^2-2x_3 x-2y_3 y-2z_3 z&=d_3^2-x_3^2-y_3^2-z_3^2\\ \end{align} however, I don't currently see this helping.

Second edit: a test case implementing Aretino's extended answer.

I must be missing a normalisation condition, can you see where I've miscalculated? One of the final distances is off...

Let $S_i(\vec m;r_i)$ be spheres with center points $m(x_i,y_i,z_i)$ and radius $r_i$.

Our known spheres are $S_1(-0.7, -0.28, 0.0; 0.7), S_2(0.5, -0.28, 0.0; 0.5), S_3(0.06, 0.45, 0.0; 0.35)$. All have no $z$ component to simplify things: this should set $\vec t = (0,0,1)$. Our unknown tangential sphere $S_0(x,y,z; 0.35)$ has a radius of $0.35$.

From these values we identify distances $d_1 = 0.7+0.35, d_2 = 0.5+0.35, d_3 = 0.35+0.35$ and solve all components to identify $\vec r=\alpha\vec u+\beta\vec v+\gamma\vec t$. Values from this example I obtain are (rounded): \begin{align} \alpha &= 0.032, &\vec u &= [-1.0, 0.0, 0.0]\\ \beta &= -0.13, &\vec v &= [-0.72, -0.69, 0.0]\\ \gamma_+ &= 0.63, &\vec t &= [0.0, 0.0, 1.0]\\ a &= -0.06, &b &= -0.10\\ c &= 0.53, &d &= -0.39\\ \vec w &= [1.4, 0.56, -0.0]\\ \end{align}

This set of values identifies an $S_0$ position that agrees with $d_1$ and $d_2$, but $d_3^{expected} = 0.7, d_3^{calculated} = 0.72$.

For completeness: the issue I had in the above example was a copy-paste error in my code. The solution provided by Aretino is correct.

• The question says that $S_1,S_2,S_3$ are known, and that the only unknown is $(x_3,y_3,z_3)$. Is that not $S_3$? Did you mean $S_0=(x_0,y_0,z_0)$ is unknown? – robjohn Jan 15 '18 at 23:48
• Yes, indeed. I've fixed that now to avoid confusion. – Geodesic Jan 16 '18 at 9:15

If I understand your question correctly, you know three vertices $S_1$, $S_2$, $S_3$ of a tetrahedron, as well as their distances $d_1$, $d_2$, $d_3$ from the fourth vertex $S_0$, whose coordinates $(x,y,z)$ must be found.

Point $S_0$ is at the intersection of three spheres, centered at $S_1$, $S_2$, $S_3$ with radii $d_1$, $d_2$, $d_3$. Its coordinates satisfy then the equations \begin{align} (x-x_1)^2+(y-y_1)^2+(z-z_1)^2=d_1^2\\ (x-x_2)^2+(y-y_2)^2+(z-z_2)^2=d_2^2\\ (x-x_3)^2+(y-y_3)^2+(z-z_3)^2=d_3^2\\ \end{align} To solve this system, subtract one of the equations from the other two, to obtain the (linear) equations of two planes. The above system is then equivalent to the quadratic system including the equations of the two planes and of one among the spheres, which should not be too difficult to solve.

That system has of course two solutions, symmetric with respect to plane $S_1S_2S_3$.

EDIT.

In order to find closed expressions of manageable size I can suggest the following approach. The system obtained by subtracting the first equation from the other two, and by the first equation itself, can be written as follows: \begin{align} \vec u\cdot\vec r=a\\ \vec v\cdot\vec r=b\\ \vec r\cdot\vec r+\vec w\cdot\vec r=c\\ \end{align} where $\vec r=(x,y,z)$ is the unknown vector, while the given quantities are: $$\vec u={(x_1-x_2,y_1-y_2,z_1-z_2)\over \sqrt{(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2}}, \quad \vec v={(x_1-x_3,y_1-y_3,z_1-z_3)\over \sqrt{(x_1-x_3)^2+(y_1-y_3)^2+(z_1-z_3)^2}},$$

$$a={d_2^2-d_1^2+x_1^2-x_2^2+y_1^2-y_2^2+z_1^2-z_2^2\over 2\sqrt{(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2}}, \quad b={d_3^2-d_1^2+x_1^2-x_3^2+y_1^2-y_3^2+z_1^2-z_3^2\over 2\sqrt{(x_1-x_3)^2+(y_1-y_3)^2+(z_1-z_3)^2}},$$

$$\vec w=(-2 x_1,-2 y_1,-2 z_1), \quad c=d_1^2-x_1^2-y_1^2-z_1^2.$$

Notice that $\vec u$ and $\vec v$ are unit vectors, and define:

$$\vec t={\vec u\times\vec v\over\sqrt{(\vec u\times\vec v)\cdot(\vec u\times\vec v)}}.$$

These vectors form a basis, so we may set

$$\vec r=\alpha\vec u+\beta\vec v+\gamma\vec t$$

and the problem of finding $\vec r$ is equivalent to that of finding $\alpha$, $\beta$ and $\gamma$. Plugging this expression for $\vec r$ into the first two equations above immediately yields: $$\alpha={a-b(\vec u\cdot\vec v)\over1-(\vec u\cdot\vec v)^2}, \quad \beta={b-a(\vec u\cdot\vec v)\over1-(\vec u\cdot\vec v)^2},$$ while the third equation becomes: $$\gamma^2+(\vec w\cdot\vec t)\gamma+d=0, \quad\hbox{where:}\quad d=\alpha^2+\beta^2+2\alpha\beta(\vec u\cdot\vec v)+ \alpha(\vec u\cdot\vec w)+\beta(\vec v\cdot\vec w)-c.$$ This last equation can finally be solved for $\gamma$: $$\gamma={1\over2}\left(-(\vec w\cdot\vec t)\pm \sqrt{(\vec w\cdot\vec t)^2-4d}\right).$$

• Thanks! This sounds like the way I'd like to go with this. Are you able to elaborate a bit more on your comment concerning the quadratic system portion of your answer? I've updated my question to show you where I've hit the wall based on my understanding of your solution. – Geodesic Jan 15 '18 at 15:32
• Looks like your goal is that of obtaining a closed form expression for $(x,y,z)$ as functions of the other variables. If so, be prepared to get very long and ugly expressions: there is no way out. If you really need them, I can feed those equations to Mathematica and post here the results. – Intelligenti pauca Jan 15 '18 at 17:07
• Ah. OK. Then that's that then. My Mathematica-foo isn't so great, so I thought when the solution in $x$ turned out to be over 400 massive terms long I had done something wrong. Any tips on an approximation method for this perhaps? – Geodesic Jan 15 '18 at 17:18
• I don't know if my edit can be of help: the idea is that of avoiding large expressions by defining suitable auxiliary quantities. – Intelligenti pauca Jan 15 '18 at 18:21
• I can confirm that the $x,y$ portions of $\vec{r}$ are correct for a number of my test cases, but seem to be obtaining negative values for $(\vec w\cdot\vec t)^2-4d$, yielding an imaginary $z$. Most likely my implementation -- I'll take a thorough look over it tomorrow. Thanks again for your help! – Geodesic Jan 15 '18 at 21:09