Using multiple solutions from equations in further calculations, part II

Question

First: This is a follow-up question to my question here, since I wasn't specific on what I really wanted.

Condider the following code:

ClearAll[a1, b1, a2, b2, xC1, yC1, xC2, yC2]
a1 := 2
b1 := 1
a2 := 1
b2 := 1/2
xC1 := 0
yC1 := 1
xC2 := 0
yC2 := 4

eqn1 = a1 (a1^2 b2^2 m^2 + b1^2 (b2^2 - (c + m xC1 - yC1)^2)) == 
       b2 (a1^2 (a2^2 m^2 + b2^2) - a2^2 (c + m xC2 - yC2)^2);
eqn2 = (a1^2 m (c - yC1) - b1^2 xC1)^2 ==
       (a1^2 m^2 + b1^2) (b1^2 xC1^2 + a1^2 (c - yC1)^2 - a1^2 b2^2);
solCM = Solve[eqn1 && eqn2, {c, m}];
Clear[x]
Evaluate[Array[x, 4]] = {(a1^2 m (c - yC1) - b1^2 xC1)/(a1^2 m^2 + b1^2)} /. solCM // Transpose // Flatten;
Simplify[x /@ Range[4]]
Clear[x]
Evaluate[Array[x, 4]] = {(a2^2 m (c - yC2) - b2^2 xC2)/(a2^2 m^2 + b2^2)} /. solCM // Transpose // Flatten;
Simplify[x /@ Range[4]]

This gives $2 \cdot 4 = 8$ values for $x$ as output, as it should. However, it is not really what I want!

I would like to insert the four $x$-values coming from

(a1^2 m (c - yC1) - b1^2 xC1)/(a1^2 m^2 + b1^2)

in the equation

((x - xC1)/a1)^2 + ((y - yC1)/b1)^2 == 1

and solve for $y$, and then give the output as four coordinates on the form $(x_i,y_i)$.

Furthermore, I would like to do a similar thing for

(a2^2 m (c - yC2) - b2^2 xC2)/(a2^2 m^2 + b2^2)

with the output values coming from the expression being inserted in

((x - xC2)/a2)^2 + ((y - yC2)/b2)^2 == 1

and solved for $y$, and give the output as another four pairs of coordinates.

P.S. The eight coordinate sets are the points at which the inner and outer tangents to two ellipses touch them. (The parameters at the beginning of the code are the major and minor axes and the coorinate sets for the centres of the two ellipses.)

Answer 1

EDIT: Comment by OP changed equation in solXY[2]. Also, removed [[1]] from Solve in SolXY to get all solutions.

ClearAll[a1, b1, a2, b2, x, y, xC1, yC1, xC2, yC2, eqn, solXY, pts]
a1 = 2;
b1 = 1;
a2 = 1;
b2 = 1/2;
xC1 = 0;
yC1 = 1;
xC2 = 0;
yC2 = 4;

eqn[1] = a1 (a1^2 b2^2 m^2 + b1^2 (b2^2 - (c + m xC1 - yC1)^2)) == 
   b2 (a1^2 (a2^2 m^2 + b2^2) - a2^2 (c + m xC2 - yC2)^2);

eqn[2] = (a1^2 m (c - yC1) - b1^2 xC1)^2 == (a1^2 m^2 + b1^2) (b1^2 xC1^2 + 
      a1^2 (c - yC1)^2 - a1^2 b2^2);

solCM = Solve[eqn[1] && eqn[2], {c, m}];

Clear[x, y]

solXY[1] = Solve[{
        x == (a1^2 m (c - yC1) - b1^2 xC1)/(a1^2 m^2 + b1^2),
        ((x - xC1)/a1)^2 + ((y - yC1)/b1)^2 == 1},
       {x, y}] /. solCM // Simplify // Flatten[#, 1] & // SortBy[#, N] &;

pts[1] = {x, y} /. solXY[1]

(*  {{-(Sqrt[35]/6), 1 - Sqrt[109]/12}, {-(Sqrt[35]/6), 
  1 + Sqrt[109]/12}, {-(Sqrt[3]/2), 1 - Sqrt[13]/4}, {-(Sqrt[3]/2), 
  1/4 (4 + Sqrt[13])}, {Sqrt[3]/2, 1 - Sqrt[13]/4}, {Sqrt[3]/2, 
  1/4 (4 + Sqrt[13])}, {Sqrt[35]/6, 1 - Sqrt[109]/12}, {Sqrt[35]/6, 
  1 + Sqrt[109]/12}}  *)

pts[1] // N

(*  {{-0.986013, 0.129974}, {-0.986013, 1.87003}, {-0.866025, 
  0.0986122}, {-0.866025, 1.90139}, {0.866025, 0.0986122}, {0.866025, 
  1.90139}, {0.986013, 0.129974}, {0.986013, 1.87003}}  *)

solXY[2] = Solve[{
        x == (a1^2 m (c - yC1) - b1^2 xC1)/(a1^2 m^2 + b1^2),
        ((x - xC2)/a2)^2 + ((y - yC2)/b2)^2 == 1},
       {x, y}] /. solCM // Simplify // Flatten[#, 1] & // SortBy[#, N] &;

pts[2] = {x, y} /. solXY[2]

(*  {{-(Sqrt[35]/6), 47/12}, {-(Sqrt[35]/6), 49/12}, {-(Sqrt[3]/2), 15/
  4}, {-(Sqrt[3]/2), 17/4}, {Sqrt[3]/2, 15/4}, {Sqrt[3]/2, 17/4}, {Sqrt[35]/6,
   47/12}, {Sqrt[35]/6, 49/12}}  *)

pts[2] // N

(*  {{-0.986013, 3.91667}, {-0.986013, 4.08333}, {-0.866025, 3.75}, {-0.866025, 
  4.25}, {0.866025, 3.75}, {0.866025, 4.25}, {0.986013, 3.91667}, {0.986013, 
  4.08333}}  *)

Verifying,

And @@ Apply[And, Table[eqn[1] && eqn[2] /. solCM /. solXY[n], {n, 2}], 2]

(*  True  *)