Common Tangent to two Ellipses
Here's a elegant solution, if I do say so myself. ;-)
The beauty is in the function:
d(α)= d1*cos(2α) + d2*sin(2α) - d3*cos(α) - d4*sin(α) + d5
I don't have any intuitive sense of the geometric reason for this, but it works.
Setup
#========================================
# Common Tangent to two Ellipses
#========================================
E1a = (-4, 1)
E1b = (-2,-1)
E1c = (-5,-1)
ellipse1 = Ellipse(E1a, E1b, E1c)
E2a = ( 3, 1)
E2b = ( 2,-1)
E2c = ( 4, 0)
ellipse2 = Ellipse(E2a, E2b, E2c)
# Transform the coordinate system (without loss of generality)
# so that ellipse1 becomes a standard circle (x^2 + y^2 = 1)
#----------------------------------------------------------------------
Cen1 = Center(ellipse1)
Maj1 = MajorAxis(ellipse1)
rot = atan2d(x(Maj1), y(Maj1)) + 180°
a1 = SemiMajorAxisLength(ellipse1)
b1 = SemiMinorAxisLength(ellipse1)
Mt = {{1,0,-x(Cen1)}, {0,1,-y(Cen1)}, {0,0,1}}
Mr = {{cos(rot),-sin(rot),0}, {sin(rot),cos(rot),0}, {0,0,1}}
Ms = {{1/a1,0,0}, {0,1/b1,0}, {0,0,1}}
Me = Ms Mr Mt
# Create matrix representation of ellipse2 (line-conic)
#----------------------------------------------------------------------
e2 = ApplyMatrix(Me, ellipse2)
c = Coefficients(e2)
M2 = {{c(1),c(4)/2,c(5)/2},{c(4)/2,c(2),c(6)/2},{c(5)/2,c(6)/2,c(3)}}
iM2= Invert(M2)
# Make a function for "signed distance" between tangent and e2
# If distance = 0 then tangent of circle is tangent of ellipse e2
#----------------------------------------------------------------------
d1 = iM2(1,1) - iM2(2,2)
d2 = 2*iM2(1,2)
d3 = 4*iM2(1,3)
d4 = 4*iM2(2,3)
d5 = iM2(1,1) + iM2(2,2) + 2*iM2(3,3)
d(α)= d1*cos(2α) + d2*sin(2α) - d3*cos(α) - d4*sin(α) + d5
Sn = Solutions(d = 0)
# Transform tangents back to the original coordinate system
#----------------------------------------------------------------------
iMe = Invert(Me)
t1 = ApplyMatrix(iMe,(cos(Sn(1))*x + sin(Sn(1))*y = 1))
t2 = ApplyMatrix(iMe,(cos(Sn(2))*x + sin(Sn(2))*y = 1))
t3 = ApplyMatrix(iMe,(cos(Sn(3))*x + sin(Sn(3))*y = 1))
t4 = ApplyMatrix(iMe,(cos(Sn(4))*x + sin(Sn(4))*y = 1))
# Some setting for created objects
#----------------------------------------------------------------------
Lobj = {"E1a", "E1b", "E1c", "ellipse1"}
Execute(Zip("SetDynamicColor("+obj+",1,0,0)", obj,Lobj))
Execute(Zip("ShowLabel("+obj+",false)", obj,Lobj))
Lobj = {"E2a", "E2b", "E2c", "ellipse2"}
Execute(Zip("SetDynamicColor("+obj+",0,0,1)", obj,Lobj))
Execute(Zip("ShowLabel("+obj+",false)", obj,Lobj))
Lobj = {"Cen1", "Maj1", "e2", "d"}
Execute(Zip("SetVisibleInView("+obj+",1,false)", obj,Lobj))
Lobj = {"t1", "t2", "t3", "t4"}
Execute(Zip("ShowLabel("+obj+",false)", obj,Lobj))
Delete(Lobj)