0
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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  *)
$\endgroup$
  • $\begingroup$ 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"? $\endgroup$ – Svend Tveskæg Apr 17 '17 at 19:24
  • $\begingroup$ @SvendTveskæg - So what have you tried? If there was a problem, what was it? $\endgroup$ – Bob Hanlon Apr 17 '17 at 19:27
  • $\begingroup$ I'm really a novice at Mathematica, so I don't know where to start / what to try. $\endgroup$ – Svend Tveskæg Apr 17 '17 at 19:30
  • $\begingroup$ 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. $\endgroup$ – Svend Tveskæg Apr 18 '17 at 18:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.