How can I reduce my solve (the number of the planes tangent to the three spheres)?

Question

Let's be three spheres $(S_1)$ has centre $A(1,2,1)$, radius $R=2$, $(S_2)$ has centre $B(3, -1, 1)$, radius $R=1$ and $(S_3)$ has centre $C(-1, -1, 1)$, radius $R=1$. How many planes tangent to the three spheres?

I suppose the equation of the plane has the form $a x + by + c z + d=0$, where $a^2 + b^2 + c^2 = 1.$

The distance from the point $A$ to the plane is $|a + 2b + c+d|$. The plane tangent to the three sphere $(S_1)$ if and only if $|a + 2b + c+d|=2$. Similarly to the spheres $(S_2)$ and $(S_3)$. My code

f[x_, y_, z_] := Abs[a x + b y + c z + d];
Solve[{f[1, 2, 1] == 2, f[3, -1, 1] == 1, f[-1, -1, 1] == 1, 
  a^2 + b^2 + c^2 == 1}, {a, b, c, d}, Reals]

I got

{{a -> -(1/2), b -> -(2/3), c -> -(Sqrt[11]/6), d -> 1/6 (-1 + Sqrt[11])}, {a -> -(1/2), b -> -(2/3), c -> Sqrt[11]/6, d -> 1/6 (-1 - Sqrt[11])}, {a -> -(1/2), b -> 2/3, c -> -(Sqrt[11]/6), d -> 1/6 (7 + Sqrt[11])}, {a -> -(1/2), b -> 2/3, c -> Sqrt[11]/6, d -> 1/6 (7 - Sqrt[11])}, {a -> 0, b -> -1, c -> 0, d -> 0}, {a -> 0, b -> -(1/3), c -> -((2 Sqrt[2])/3), d -> 2/3 (-2 + Sqrt[2])}, {a -> 0, b -> -(1/3), c -> (2 Sqrt[2])/3, d -> -(2/3) (2 + Sqrt[2])}, {a -> 0, b -> 1/3, c -> -((2 Sqrt[2])/3), d -> 2/3 (2 + Sqrt[2])}, {a -> 0, b -> 1/3, c -> (2 Sqrt[2])/3, d -> -(2/3) (-2 + Sqrt[2])}, {a -> 0, b -> 1, c -> 0, d -> 0}, {a -> 1/2, b -> -(2/3), c -> -(Sqrt[11]/6), d -> 1/6 (-7 + Sqrt[11])}, {a -> 1/2, b -> -(2/3), c -> Sqrt[11]/6, d -> 1/6 (-7 - Sqrt[11])}, {a -> 1/2, b -> 2/3, c -> -(Sqrt[11]/6), d -> 1/6 (1 + Sqrt[11])}, {a -> 1/2, b -> 2/3, c -> Sqrt[11]/6, d -> 1/6 (1 - Sqrt[11])}

Now I count by hand and DeleteCases.

Answer 1

I think the equations are not quite what is wanted.

For each sphere_j there is a point {x[j],y[j],z[j]} that must satisfy the equations that (i) it is on that particular sphere (ii) it is on the plane (iii) the gradient of the sphere, at that point, is parallel to the plane direction vector {a,b,c}. We also use a normalization equation for the plane parameters. The code below will set up this system.

centers = {{1, 2, 1}, {3, -1, 1}, {-1, -1, 1}};
radii = {2, 1, 1};
sphereVals = Thread[{centers, radii}];

plane[x_, y_, z_] := {a, b, c, d}.{x, y, z, -1}
spheres = 
  MapIndexed[
   With[{coords = {x[#2[[1]]], y[#2[[1]]], z[#2[[1]]]} - #[[1]]}, 
     coords.coords - #[[2]]^2] &, sphereVals];
grads = MapIndexed[Grad[#, {x[#2[[1]]], y[#2[[1]]], z[#2[[1]]]}] &, 
   spheres];
planeHitsSpherePolys = Map[plane[x[#], y[#], z[#]] &, Range[3]];
planeTangSpherePolys = Flatten[MapIndexed[#1 - k[#2[[1]]]*{a, b, c} &, grads]];

So here are the equations.

allPolys = 
 Join[spheres, planeHitsSpherePolys, 
  planeTangSpherePolys, {a^2 + b^2 + c^2 + d^2 - 1}]

(* {-4 + (-1 + x[1])^2 + (-2 + y[1])^2 + (-1 + 
    z[1])^2, -1 + (-3 + x[2])^2 + (1 + y[2])^2 + (-1 + 
    z[2])^2, -1 + (1 + x[3])^2 + (1 + y[3])^2 + (-1 + z[3])^2, -d + 
  a x[1] + b y[1] + c z[1], -d + a x[2] + b y[2] + c z[2], -d + 
  a x[3] + b y[3] + c z[3], -a k[1] + 2 (-1 + x[1]), -b k[1] + 
  2 (-2 + y[1]), -c k[1] + 2 (-1 + z[1]), -a k[2] + 
  2 (-3 + x[2]), -b k[2] + 2 (1 + y[2]), -c k[2] + 
  2 (-1 + z[2]), -a k[3] + 2 (1 + x[3]), -b k[3] + 
  2 (1 + y[3]), -c k[3] + 2 (-1 + z[3]), -1 + a^2 + b^2 + c^2 + d^2} *)

We don't really care about the ratios between direction vector and gradients, and from the wording of the question I will guess we also do not require the three intersection points.

solns = {a, b, c, d} /. Solve[allPolys == 0];

Numerically:

N[solns]

    (* {{0., -1., 0., 0.}, {0., 1., 0., 0.}, {0., 
  0.134077, -0.379228, -0.915538}, {0., -0.134077, 0.379228, 
  0.915538}, {0., -0.310496, -0.878216, 0.363769}, {0., 0.310496, 
  0.878216, -0.363769}, {-0.251372, 
  0.335162, -0.277902, -0.864436}, {0.251372, -0.335162, 0.277902, 
  0.864436}, {0.405875, 
  0.541167, -0.448712, -0.584004}, {-0.405875, -0.541167, 0.448712, 
  0.584004}, {0.426112, -0.568149, -0.471085, 0.523177}, {-0.426112, 
  0.568149, 
  0.471085, -0.523177}, {-0.46644, -0.62192, -0.515669, -0.360189}, \
{0.46644, 0.62192, 0.515669, 0.360189}} *)

These come in pairs since the negative of a solution is also a solution (giving the same plane of course). So taking that into account brings it to seven such planes, as was claimed in the original post. Also we can renormalize and simplify. The rationale for renormalizing is that maybe the original norm caused an excess of radicals, and dividing each solution by the same nonzero component could give a better form.

possolns = Select[solns, #[[2]] >= 0 &];
simpler = Simplify[Map[#/Max[Abs[#]] &, possolns]]

(* Out[79]= {{0, 1, 0, 0}, {0, 1/(4 + 2 Sqrt[2]), 1 - Sqrt[2], -1}, {0, 
  1/(2 Sqrt[2]), 1, 1 - Sqrt[2]}, {-(3/(7 + Sqrt[11])), 4/(
  7 + Sqrt[11]), 1/38 (11 - 7 Sqrt[11]), -1}, {3/(1 + Sqrt[11]), 4/(
  1 + Sqrt[11]), 1/10 (-11 + Sqrt[11]), -1}, {-(3/4), 1, Sqrt[11]/4, 
  1/4 (-7 + Sqrt[11])}, {3/4, 1, Sqrt[11]/4, 1/4 (-1 + Sqrt[11])}} *)

Someone with better graphics chops than myself might be able to show all these planes along with the spheres.