# Using multiple solutions from equations in further calculations, part II

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.)

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  *)

• Hmmm! I've missed something; I can see that all eight point have to calculated from only (a1^2 m (c - yC1) - b1^2 xC1)/(a1^2 m^2 + b1^2) (and not half of them from (a2^2 m (c - yC2) - b2^2 xC2)/(a2^2 m^2 + b2^2)). Can the code them be made "simpler"? Commented Apr 17, 2017 at 19:24
• @SvendTveskæg - So what have you tried? If there was a problem, what was it? Commented Apr 17, 2017 at 19:27
• I'm really a novice at Mathematica, so I don't know where to start / what to try. Commented Apr 17, 2017 at 19:30
• I was mistaken (again); it really has to be both (a1^2 m (c - yC1) - b1^2 xC1)/(a1^2 m^2 + b1^2) and (a2^2 m (c - yC2) - b2^2 xC2)/(a2^2 m^2 + b2^2). Nevermind, since I've managed to put together all my different pieces of code in order to make the "program" work. Commented Apr 18, 2017 at 18:28