Circle tangent to ellipse

Question

A small circle is tangent to the lower right quarter of a larger ellipse. The ellipse is straight (not rotated), and its lowest point is centered at origo.

The following is given: "a", the axis of the ellipse along the x axis. "b", the axis of the ellipse along the y axis. "r", the radius of the circle. "h", the top y coordinate of the circle.

I am seeking the x coordinate of the tangent point between the ellipse and the circle as a function of a, b, r and h.

Answer 1

The ellipse has the equation

$ \dfrac{x^2}{a^2} + \dfrac{(y - b)^2}{b^2} = 1 \hspace{20pt}(1)$

Now let the circle have a center $(p, q)$, then $q = h - r$, and its equation is

$ (x - p)^2 + (y - q)^2 = r^2 \hspace{20pt}(2)$

The unknowns here are $x, y, p$

In addition to the above two equations, we know that the normal to the ellipse at $(x,y)$ and the normal to the circle at the same point, are aligned.

The normal vector (also called gradient) of a function of two variables, $F(x,y) = 0$ is given by $N = ( \dfrac{\partial F}{\partial x} , \dfrac{\partial F}{\partial y})$.

Hence, the normal vector to the ellipse is

$N_1 = (N_{1x}, N_{1y}) = \left( \dfrac{2 x}{ a^2} , \dfrac{2 (y - b)}{b^2} \right) $

similarly, the normal vector to the circle is

$N_2 = (N_{2x}, N_{2y}) = ( 2 (x - p) , 2 (y - q) ) $

If $N_1$ is aligned with $N_2$ then

$ N_1 = k N_2 $

where $k$ is a constant.

Hence,

$(N_{1x} , N_{1y} ) = k (N_{2x}, N_{2y} ) $

From which,

$ k = \dfrac{N_{1x}}{N_{2x}} = \dfrac{N_{1y}}{N_{2y}} $

Using cross multiplication on the equality on the right of the above equation, we obtain,

$ N_{1x} N_{2y} = N_{2x} N_{1y} $

Substituting the expressions for $N_{1x}, N_{1y}, N_{2x}, N_{2y}$, gives us,

$ \dfrac{2 x}{ a^2} \cdot 2 (y - q) = 2 (x - p) \cdot \dfrac{2 (y - b)}{b^2} $

Dividing the whole equation by $4$ to cancel the $2$'s, and then multiplying through by $a^2 b^2$ gives us,

$ b^2 x (y - q) = a^2 (x - p) (y - b) \hspace{20pt}(3)$

Now, we'll to eliminate $x$ and $(x - p)$ from equation $(3)$

First square equation $(3)$, this gives us,

$ b^4 x^2 (y - q)^2 = a^4 (x - p)^2 (y - b)^2 \hspace{20pt}(4) $

From equation $(1)$,

$ x^2 = a^2 \left(1 - \dfrac{(y - b)^2}{b^2} \right) \hspace{20pt}(5) $

And, from equation $(2)$,

$ (x - p)^2 = r^2 - (y - q)^2 \hspace{20pt}(6)$

Substitute $(5)$ and $(6)$ into $(4)$. This gives

$ b^4 (y - q)^2 a^2 \left(1 - \dfrac{(y - b)^2}{b^2} \right) = a^4 (y - b)^2 (r^2 - (y - q)^2 ) \hspace{20pt} (7) $

This simplifies to

$ b^2 (y - q)^2 ( - y^2 + 2 y b ) = a^2 (y - b)^2 (r^2 - (y - q)^2 ) \hspace{20pt}(8)$

Expanding equation $(8)$, yields the quartic polynomial equation:

$ a_4 y^4 + a_3 y^3 + a_2 y^2 + a_1 y + a_0 = 0 \hspace{20pt} (9)$

where

$a_4 = b ^ 2 - a ^ 2 $

$a_3 = 2 (b + q) (a ^ 2 - b ^ 2) $

$a_2 = a ^ 2 (r ^ 2 - q ^ 2 - b ^ 2 - 4 b q) - b ^ 2 (-q ^ 2 - 4 b q) $

$a_1 = a ^ 2 (2 b (q ^ 2 - r ^ 2) + 2 q b ^ 2) - 2 b ^ 3 q ^ 2$

$a_0 = a ^ 2 b ^ 2 (r ^ 2 - q ^ 2) $

The roots of this quartic polynomial can be solved by standard routines in any math library. These routines find the exact value of the root using exact formulas for the quartic polynomial roots that depend only on the coeffcients $a_4$ through $a_0$. The exact way these roots are calculated is beyond the scope of this answer. The interested reader is referred to this Wikipedia page.

To select the correct root, note that

$ x = a \sqrt{ 1 - \dfrac{(y - b)^2}{b^2} } \hspace{20pt}(10) $

Because we know that $ x \gt 0 $, while from equation $(6)$, and since $p \gt x$,

$ (x - p) = - \sqrt{ r^2 - (y - q)^2 } \lt 0 \hspace{20pt}(6)$

Now look at equation $(3)$

$ b^2 x (y - q) = a^2 (x - p) (y - b) $

Since we want tangency below the line $y = b $, then we firstly require that $ (y - b) \lt 0$, this will make the right hand side positive. Now looking at the left hand side, we see that $x \gt 0$, so we must require that $ y \gt q $. Hence, in summary we will select $y$ such that $ q \lt y \lt b $. Once we have found $y$, then the calculation of $x$ and $p$ follows as stated above.

Here is an example, where $a = 4, b = 8, r = 2 , h = 5 $.

Running a small code which calls a quartic polynomial root finder, and after selecting the appropriate $y$, the program gives the tangent point as $(x,y)$ and center of circle $(p, q)$ with

$ x = 3.34773271274127 $

$ y = 3.62155932595085 $

$ p = 5.24869667987045 $

$ q = 3 $

Answer 2

The ellipse $E_{x,y}\to \left(\frac xa\right)^2+\left(\frac{y-b}{b}\right)^2=1$ in polar coordinates reads

$$ E_{\theta}\to\cases{x=a\cos\theta\\ y = b+b\sin\theta} $$

The normal to $E_{x,y}$ is given by

$$ \vec n_{x,y}= \left(\frac{x}{a^2 \sqrt{\frac{x^2}{a^4}+\frac{(b-y)^2}{b^4}}},\frac{y-b}{b^2 \sqrt{\frac{x^2}{a^4}+\frac{(b-y)^2}{b^4}}}\right) $$

and consequently

$$ \vec n_{\theta} = \left(\frac{\cos \theta }{a \sqrt{\frac{\cos ^2\theta }{a^2}+\frac{\sin ^2\theta }{b^2}}},\frac{\sin \theta }{b \sqrt{\frac{\cos ^2\theta }{a^2}+\frac{\sin ^2\theta }{b^2}}}\right) $$

then

$$ \mathcal{E}_{\theta} = \left(a \cos\theta,b+b\sin\theta\right)+ r \vec n_{\theta} $$

describes a set of points expanded from $E_{\theta}$ in the normal direction, with length $r$.

Now $\mathcal{E}_{\theta}$ intersects with the line $y = h - r$ where the circle center should be located, or better:

$$ b+b\sin\theta+\frac{r\sin \theta }{b\sqrt{\frac{\cos ^2\theta }{a^2}+\frac{\sin ^2\theta }{b^2}}}= h-r $$

or

$$ f(\theta)=b(1+\sin\theta)+\frac{r a\sin\theta}{\sqrt{b^2-(b^2-a^2)\sin^2\theta}}-h+r=0 $$

now choosing an iterative process like Newton's to approximate an initial guess,

$$ \theta_{k+1} = \theta_k - f(\theta_k)/f'(\theta_k) $$

we are done.

Taking $a=1,b=2,h=1/3, r=1/5$ and $\theta_0=0$ in three iterations we obtain $\theta = -0.8675956$. Follows a graphic showing $E_{\theta}$ yellow continuous, $\mathcal{E}_{\theta}$ yellow dashed, $y=h-r$ in blue dashed and in red the circle center calculated as $(0.647982738, 0.133333333)$.

Attached a python script to perform the calculations

import math

a = 1
b = 2
h = 1./3.
r = 0.2
theta0 = 0.
n = 7
error = 0.0000001

def Phi(theta):
    s = math.sin(theta)
    c = math.cos(theta)
    p1 = b**2-(b**2 - a**2)*s**2
    p2 = a*r/math.sqrt(p1) + b
    p3 = s*p2 + b - h + r
    p4 = c*p2 - a*r*(a**2 - b**2)*s**2*c/math.sqrt(p1)**3
    return theta - p3/p4


for k in range(n):
    theta1 = Phi(theta0)
    if abs(theta1-theta0) < error:
        print("result in ", k, "iterations. theta  = ", theta1)
        c = math.cos(theta1)
        s = math.sin(theta1)
        x0 = a*c + r*c*b/math.sqrt(a**2*s**2 + b**2*c**2)
        y0 = b + s*(b + a*r/math.sqrt(a**2*s**2 + b**2*c**2))
        print("circle center = (", x0, y0, ")")
        break
    theta0 = theta1