Given an ellipse and a reference point, how to find the two lines that are tangent to the ellipse?

Question

I have an ellipse, possibly rotated and shifted from the origin, which is given by a parametrization similar to this one: $$ \begin{aligned} x &= x_0 + a\cos\theta\cos\alpha - b\sin\theta\sin\alpha \\ y &= y_0 + a\cos\theta\sin\alpha + b\sin\theta\cos\alpha \\ \end{aligned} $$

where $(x_0,y_0)$ is the center of the ellipse, $(a,b)$ are its semi-axes, $\alpha$ is the rotation angle and $\theta\in [0,2\pi]$.

The reference point is given by $(p_0,q_0)$. Now I want to determine the parameters of the lines $y = mx + h$ that are tangent to the ellipse and go through that reference point. To do so, I compute $h$ from the reference point, $h = q_0-mp_0$, and then plug the ellipse equations into the line equation: $$ \underbrace{y_0-q_0}_{\equiv q} + a\cos\theta\sin\alpha + b\sin\theta\cos\alpha = m\left(\underbrace{x_0-p_0}_{\equiv p} + a\cos\theta\cos\alpha - b\sin\theta\sin\alpha \right) $$ From this I can get an equation for the slope $m$ in dependence on $\theta$. I also know that the slope of a tangent to the ellipse is given by $\frac{y'(\theta)}{x'(\theta)}$, i.e.: $$ m = \frac{a\sin\theta\sin\alpha - b\cos\theta\cos\alpha}{a\sin\theta\cos\alpha + b\cos\theta\sin\alpha} $$

Now, if I equate the two above equations for $m$, I get an expression which can be solved for $\theta$. I would expect that there are two possible solutions to this equations, as there are two tangent lines. I used WolframAlpha to solve this equation (see here), but it gave me four different solutions (which differ in the combination of two signs); I'm including a screenshot of one of the solutions below since it's quite long ($x\equiv\theta, t\equiv\alpha$):

The four solutions emerge as $x = \pm\cos^{-1}((\pm\sqrt{\dots}$. There are no additional constraints mentioned, so something must have gone wrong, but I can't figure out what. I verified that two of those four solutions are valid with the help of some examples (testing various reference points p0,q0):

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import numpy as np
from numpy import arccos, cos, pi, sign, sin, sqrt

a, b = 3, 1  # ellipse parameters
center = 1, 2  # center of ellipse
angle = 60  # rotation of ellipse (counter-clockwise) [deg]
p0, q0 = -2, 3  # the reference point (x-,y-coordinates)

p = center[0] - p0
q = center[1] - q0
t = angle/180*pi
thetas = [
    s1*arccos((s2*sqrt(a**2*b**4*(2*p*cos(t)**3 + 2*p*sin(t)**2*cos(t) + 2*q*sin(t)**3 + 2*q*sin(t)*cos(t)**2)**2 + 4*a**2*(-b**2*sin(t)**4 - b**2*cos(t)**4 - 2*b**2*sin(t)**2*cos(t)**2 + p**2*sin(t)**2 - 2*p*q*sin(t)*cos(t) + q**2*cos(t)**2)*(a**2*p**2*sin(t)**2 - 2*a**2*p*q*sin(t)*cos(t) + a**2*q**2*cos(t)**2 + b**2*p**2*cos(t)**2 + 2*b**2*p*q*sin(t)*cos(t) + b**2*q**2*sin(t)**2)) - a*b**2*(2*p*cos(t)**3 + 2*p*sin(t)**2*cos(t) + 2*q*sin(t)**3 + 2*q*sin(t)*cos(t)**2))/(2*(a**2*p**2*sin(t)**2 - 2*a**2*p*q*sin(t)*cos(t) + a**2*q**2*cos(t)**2 + b**2*p**2*cos(t)**2 + 2*b**2*p*q*sin(t)*cos(t) + b**2*q**2*sin(t)**2)))
    for s1 in [-1, 1] for s2 in [-1, 1]
]

fig, ax = plt.subplots()
ax.add_patch(Ellipse(center, 2*a, 2*b, angle, color='tab:blue', alpha=0.4))
ax.scatter([p0], [q0], color='tab:blue', s=50)
for theta in thetas:
    ex = center[0] + a*cos(theta)*cos(t) - b*sin(theta)*sin(t)
    ey = center[1] + a*cos(theta)*sin(t) + b*sin(theta)*cos(t)
    ax.axline((p0,q0), (ex,ey), color='tab:orange', lw=1)
    ax.scatter([ex], [ey], color='tab:orange', s=25)
plt.show()

So it seems that all that's left is to find a way how to pick the right combinations of signs for the two tangent points. However, I don't know what condition I should apply. I tried to discriminate based on the quadrant of $(p_0-x_0,q_0-y_0)$, but that didn't work. Also, since there are 6 different ways to pick 2 out of 4 possible sign combinations, it seems that the condition should include 6 different "sections/intervals".

Answer 1

One way is to use the equation which is $$(\frac{(x-x_0)\cos(\alpha)+(y-y_0)\sin(\alpha)}{a})^2+(\frac{-(x-x_0)\sin(\alpha)+(y-y_0)\cos(\alpha)}{b})^2-1=0.$$

Then the equation of the line pair (tangents from $(p_0,q_0)$) are by joachimsthal $s\cdot s_{11}-s_1^2=0$ or

$((\frac{(x-x_0)\cos(\alpha)+(y-y_0)\sin(\alpha)}{a})^2+(\frac{-(x-x_0)\sin(\alpha)+(y-y_0)\cos(\alpha)}{b})^2-1)((\frac{(p_0-x_0)\cos(\alpha)+(q_0-y_0)\sin(\alpha)}{a})^2+(\frac{-(p_0-x_0)\sin(\alpha)+(q_0-y_0)\cos(\alpha)}{b})^2-1)\\-((\sin(\alpha)^2/b^2+\cos(\alpha)^2/a^2)xp_0+((\cos(\alpha)\sin(\alpha))/a^2-(\cos(\alpha)\sin(\alpha))/b^2)(xq_0+p_0y)+(\cos(\alpha)^2/b^2+\sin(\alpha)^2/a^2)yq_0+((\cos(\alpha)\sin(\alpha)y_0)/b^2-(\cos(\alpha)\sin(\alpha)y_0)/a^2-(\sin(\alpha)^2x_0)/b^2-(\cos(\alpha)^2x_0)/a^2)(x+p_0)+((-(\cos(\alpha)^2y_0)/b^2)-(\sin(\alpha)^2y_0)/a^2+(\cos(\alpha)\sin(\alpha)x_0)/b^2-(\cos(\alpha)\sin(\alpha)x_0)/a^2)(y+q_0)+((\cos(\alpha)^2y_0^2)/b^2+(\sin(\alpha)^2y_0^2)/a^2-(2\cos(\alpha)\sin(\alpha)x_0y_0)/b^2+(2\cos(\alpha)\sin(\alpha)x_0y_0)/a^2+(\sin(\alpha)^2x_0^2)/b^2+(\cos(\alpha)^2x_0^2)/a^2-1))^2=0$

which simplifies to

$$\frac{((y_0-q_0)^2-((b^2-a^2)\cos(2α)+a^2+b^2)/2)(x-p_0)^2+(\sin(2α)(a^2-b^2)-2(x_0-p_0)(y_0-q_0))(x-p_0)(y-q_0)+((x_0-p_0)^2+((b^2-a^2)\cos(2α)-(a^2+b^2))/2)(y-q_0)^2}{a^2b^2}=0$$ which you can factor by setting $t=\frac{x-p_0}{y-q_0}$ (or its inverse) and using the quadratic formula to get

$$(x-p_0)-(y-q_0)\frac{2(x_0-p_0)(y_0-q_0)+\sin(2α)(b^2-a^2)\pm\sqrt{D}}{2(y_0-q_0)^2-((b^2-a^2)\cos(2α)+b^2+a^2)}=0$$

where

$$D=2(a^2+b^2+(a^2-b^2)\cos(2\alpha))(y_0-q_0)^2 +4(b^2-a^2)\sin(2\alpha)(x_0-p_0)(y_0-q_0) +2((b^2-a^2)\cos(2\alpha)+b^2+a^2)(x_0-p_0)^2-4a^2b^2$$

Answer 2

Let $x=x_0+a\cos\theta\cos\alpha-b\sin\theta\sin\alpha$, $y=y_0+a\cos\theta\sin\alpha+b\sin\theta\cos\alpha$ is ellipse point with parameter $\theta\in[0;2\pi)$. Then slope of line going through this point and reference point $(p_0,q_0)$ is $$k=\frac{y-q_0}{x-p_0}=\frac{y_0-q_0+a\cos\theta\sin\alpha+b\sin\theta\cos\alpha}{x_0-p_0+a\cos\theta\cos\alpha-b\sin\theta\sin\alpha}=\frac{A+B\cos\theta+C\sin\theta}{D+E\cos\theta+F\sin\theta}$$

$$\frac{dk}{d\theta}=\frac{CE-BF+(CD-AF)\cos\theta+(AE-BD)\sin\theta}{(D+E\cos\theta+F\sin\theta)^2}.$$

Tangent points will be determined by $$CE-BF+(CD-AF)\cos\theta+(AE-BD)\sin\theta=0$$ $$\frac{(CD-AF)\cos\theta+(AE-BD)\sin\theta}{\sqrt{(CD-AF)^2+(AE-BD)^2}}=\frac{BF-CE}{\sqrt{(CD-AF)^2+(AE-BD)^2}}$$ $$\cos\beta=\frac{AE-BD}{\sqrt{(CD-AF)^2+(AE-BD)^2}},\; \sin\beta=\frac{CD-AF}{\sqrt{(CD-AF)^2+(AE-BD)^2}}$$ $$\sin(\theta+\beta)=\frac{BF-CE}{\sqrt{(CD-AF)^2+(AE-BD)^2}}$$ $$\theta_1=-\beta+\arcsin\frac{BF-CE}{\sqrt{(CD-AF)^2+(AE-BD)^2}},$$ $$\theta_2=\pi-\beta-\arcsin\frac{BF-CE}{\sqrt{(CD-AF)^2+(AE-BD)^2}}$$

The calculation procedure: given $x_0$, $y_0$, $a$, $b$, $\alpha$, $p_0$, $q_0$. Calculate $A$, $B$, $C$, $D$, $E$, $F$. Then calculate $\cos\beta$, $\sin\beta$, $\sin(\theta+\beta)$ by formulae from above. Then find $\beta$ and two values of $\theta$. Then calculate $(x,y)$ for both tangent points, using $\theta_1$ and $\theta_2$.

Answer 3

There can be only two solutions.

You can make the problem much simpler by transforming the ellipse to a circle (rotate by $-\alpha$ and stretch by $(\frac1a,\frac1b)$), and move it to the origin. Apply the same transformation to the external point.

Now you can solve by simple trigonometry or express that the point of contact defines a right angle:

$$(p-\cos t)\cos t+(q-\sin t)\sin t=0$$

or $$\cos\left(t-\arctan\frac qp\right)=\frac1{\sqrt{p^2+q^2}}.$$