Circle tangent to rotated ellipse and horizontal line

Question

I would like to find the position for the center of a circle $(x_0, y_0)$ that is tangent to both an ellipse and a horizontal line.
The ellipse is positioned at $(0,0)$ and is defined by major axis $a$ and minor axis $b$ and can be rotated by $\alpha$.
The line $h$ is horizontal and can be moved on the $y$-axis.

The Python snippet in this answer works well but doesn't allow for the ellipse to be rotated.

Edit: I'm a technical artist in the game industry working on a procedural tool that involves calculating the center of a circle to create an almost accurate 3D model. I need to understand the math behind it in order to generate code to solve it.
My math knowledge ranges from high school level to math applicable to basic 3D graphics.

Answer 1

From

$$ \left(\frac xa\right)^2+\left(\frac yb\right)^2-1=0 $$

changing of coordinates by rotating $-\theta$ we have

$$ \frac{(X \cos (\theta )+Y \sin (\theta ))^2}{a^2}+\frac{(Y \cos (\theta )-X \sin (\theta ))^2}{b^2}-1=0 $$

or

$$ \left(\frac{\cos ^2(\theta )}{a^2}+\frac{\sin ^2(\theta )}{b^2}\right)X^2+\left(\frac{1}{a^2}-\frac{1}{b^2}\right) \sin (2 \theta )XY+\left(\frac{\sin ^2(\theta )}{a^2}+\frac{\cos ^2(\theta )}{b^2}\right)Y^2-1=0 $$

and the normal vector $\vec n$ is given by

$$ \vec n = \left\{ \begin{array}{c} \frac{ \cos ^2(\theta ) (Y \tan (\theta )+X)}{a^2}+\frac{ \sin (\theta ) (X \sin (\theta )-Y \cos (\theta ))}{b^2} \\ \frac{ \sin (\theta ) (X \cos (\theta )+Y \sin (\theta ))}{a^2}+\frac{ \cos (\theta ) (Y \cos (\theta )-X \sin (\theta ))}{b^2} \\ \end{array} \right. $$

now changing to polar coordinates $$ E(t)=\cases{ X = \rho(t)\cos t\\ Y = \rho(t)\sin t } $$

we have

$$ \cases{ \rho ^2(t) \left(\frac{\cos ^2(t-\theta )}{a^2}+\frac{\sin ^2(t-\theta )}{b^2}\right)-1=0\\ \vec n(t) = \left\{ \begin{array}{c} \frac{\rho(t) \left(\left(b^2-a^2\right) \cos (t-2 \theta )+\left(a^2+b^2\right) \cos (t)\right)}{a^2 b^2} \\ \frac{\rho(t) \left(\left(a^2+b^2\right) \sin (t)+(a-b) (a+b) \sin (t-2 \theta )\right)}{a^2 b^2} \\ \end{array} \right. } $$

now taking

$$ \rho(t) = \frac{ab\sqrt{2}}{\sqrt{\left(b^2-a^2\right) \cos (2 (t-\theta ))+a^2+b^2}} $$

and substituting into

$$ \mathcal{E}_r(t)=(\rho(t)\cos t,\rho(t)\sin t) + r\frac{\vec n}{\|\vec n\|} $$

where $\mathcal{E}_r(t)$ is $E(t)$ expanded by $r$, we get at

$$ \mathcal{E}_r(t)=\frac{ab\sqrt{2}}{\sqrt{\left(b^2-a^2\right) \cos (2 (t-\theta ))+a^2+b^2}}\left(\cos t, \sin t\right)+r\frac{\vec n}{\|\vec n\|} $$

where

$$ \vec n = \left\{ \begin{array}{c} \frac{\sqrt{2} \left(\left(b^2-a^2\right) \cos (t-2 \theta )+\left(a^2+b^2\right) \cos (t)\right)}{a b \sqrt{\left(b^2-a^2\right) \cos (2 (t-\theta ))+a^2+b^2}} \\ \frac{\sqrt{2} \left(\left(a^2+b^2\right) \sin (t)+(a-b) (a+b) \sin (t-2 \theta )\right)}{a b \sqrt{\left(b^2-a^2\right) \cos (2 (t-\theta ))+a^2+b^2}} \\ \end{array} \right. $$

now, if the horizontal line is given by $Y=h$, the tangent circle center is located at the solutions for

$$ (0,1)\cdot \mathcal{E}_r(t) = h + r $$

where $r$ is the circle radius.

NOTE

Assuming the parametric values $a=1,b=3,h=1,r=0.5,\theta = -\frac{\pi}{6}$ we have

$$ (0,1)\cdot \mathcal{E}_r(t) - (h + r) = \frac{3 \sin (t)}{\sqrt{2 \sqrt{3} \sin (2 t)+2 \cos (2 t)+5}}+\frac{3 \sin (t)+2 \sqrt{3} \cos (t)}{3 \sqrt{2 \sqrt{3} \sin (2 t)+2 \cos (2 t)+5} \sqrt{\frac{4 \left(3 \sin (t)+2 \sqrt{3} \cos (t)\right)^2}{9 \left(2 \sqrt{3} \sin (2 t)+2 \cos (2 t)+5\right)}+\frac{4 \left(2 \sqrt{3} \sin (t)+7 \cos (t)\right)^2}{9 \left(2 \sqrt{3} \sin (2 t)+2 \cos (2 t)+5\right)}}}-\frac{3}{2}=0 $$

to solve for $t$. Follows a plot shoving in blue $E(t)$ in dashed light brown $\mathcal{E}_{\frac 12}(t)$ also in dashed red $Y=1+\frac 12$. The solution for $t$ can be obtained by binary search.

Follows a python script to find the solution angles and circle centers.

import numpy as np
import matplotlib.pyplot as plt


n     = 150
x0    = np.linspace(0,np.pi,n)
err   = 0.00001
a     = 1
b     = 3
h     = 1
r     = 1/2
theta =-np.pi/6

def expanded_ellipse(t):
    st   = np.sin(t)
    ct   = np.cos(t)
    stt  = np.sin(t+theta)
    ctt  = np.cos(t+theta)
    c2tt = np.cos(2*(t+theta))
    st2t = np.sin(t+2*theta)
    ct2t = np.cos(t+2*theta)
    den1 = np.sqrt((b*ctt)**2 + (a*stt)**2)
    den2 = np.sqrt(2)*np.sqrt(a**4 + b**4 + (b**4 - a**4)*c2tt)
    numx = ((a**2 + b**2)*ct + (b**2 - a**2)*ct2t)*r
    numy = ((a**2 + b**2)*st + (a**2 - b**2)*st2t)*r
    return [a*b*ct/den1 + numx/den2, a*b*st/den1 + numy/den2]


def residual(t):
    st   = np.sin(t)
    stt  = np.sin(t+theta)
    ctt  = np.cos(t+theta)
    c2tt = np.cos(2*(t+theta))
    st2t = np.sin(t+2*theta)
    num  = r*((a**2 + b**2)*st + (a**2 - b**2)*st2t)
    den  = np.sqrt(2)*np.sqrt(a**4 + b**4 + (b**4 - a**4)*c2tt)
    return a*b*st/np.sqrt((b*ctt)**2+(a*stt)**2) + num/den - h - r


def minus_residual(t):
    return -residual(t)


def binary_search(f, domain, MAX = 1000):
    start, end = domain
    if start >= end:
        raise ValueError("Domain is empty")
    mid = (start + end) / 2
    fmid = f(mid)
    step = 0
    while abs(fmid) > err and step < MAX:
        if fmid < 0:
            start = mid
        else:
            end = mid
        mid = (start + end) / 2
        fmid = f(mid)
        step += 1
    return round(mid, 7)


y0 = []

for x in x0:
    y0.append(residual(x))

plt.xlim([0,np.pi])
plt.ylim([min(y0),max(y0)])
plt.plot(x0,y0,'-',color = 'black')
plt.plot(x0,n*[0])

t1 = binary_search(residual,[1,1.7])
t2 = binary_search(minus_residual,[2,3])

print(t1, expanded_ellipse(t1))
print(t2, expanded_ellipse(t2))

NOTE

The binary_search procedure needs that the function be increasing in the searching domain: so we have residual and minus_residual.