Circle tangent to ellipse

from sympy import *
from sympy.abc import *


def main():
    # inputs
    ai = 539.2645225464834
    bi = 5255.824762975057
    ri = 3
    hi = 20
    # ---

    qi = hi - ri

    a4 = b ** 2 - a ** 2
    a3 = 2 * (b + q) * (a ** 2 - b ** 2)
    a2 = a ** 2 * (r ** 2 - q ** 2 - b ** 2 - 4 * b * q) - b ** 2 * (-q ** 2 - 4 * b * q)
    a1 = a ** 2 * (2 * b * (q ** 2 - r ** 2) + 2 * q * b ** 2) - 2 * b ** 3 * q ** 2
    a0 = a ** 2 * b ** 2 * (r ** 2 - q ** 2)
    fy = a4 * y ** 4 + a3 * y ** 3 + a2 * y ** 2 + a1 * y + a0
    solutions = solve(fy, y, dict=True)
    # the third (index 2) element of the solutions is apparently always the one we're looking for
    ye = solutions[2][y]
    # Strip away complex part (very small)
    y1 = ye.evalf(subs={a: ai, b: bi, r: ri, q: qi}).args[0]
    xe = a * (1 - (ye - b) ** 2 / b ** 2) ** (1 / 2)
    x1 = xe.evalf(subs={a: ai, b: bi, r: ri, q: qi}).args[0]
    pe = xe + (r ** 2 - (ye - q) ** 2) ** (1 / 2)
    p1 = pe.evalf(subs={a: ai, b: bi, r: ri, q: qi}).args[0]
    print(y1, x1, p1)


if __name__ == '__main__':
    main()