How to calculate intersection point between lines and ellipses using only math equation

I want to calculate the intersection points between an ellipse and a line using python. However how do I do it using only mathematical quadratic equation, when given both end points of the line and the center point and radii of the ellipse?

None of the solutions on stackoverflow describe the solution I'm looking for. They use either numpy, scipy and other libraries or use trigonometry - none of them use a purely mathematical process.

Comparing the solution in this post for example:

You have the ellipse formula:

(x - center_x)^2 / a^2 + (y - center_y)^2 / b^2 = 1

And the line formula:

y = m*x + c

The solution in linked post describes the quadratic formula derived from it on rearrangement:

A=a2m2+b2

B=2a2m(c-k)-2b2h

C=b2h2+a2(c-k)2-a2b2

I want to solve this equation using only math, solving for x, and afterwards solving for y using the line formula.

However when using above formula, the resulting intersection points are erroneous and do not lie on the ellipse or the line.

This is the code I tried:

from math import *
def GetIntersectionPoints(line_point_a, line_point_b, center_x, center_y, rad_x, rad_y):
    # use line formula to calculate m and c:
    x1, y1, x2, y2 = *line_point_a, *line_point_b
    m = (y1 - y2) / (x1 - x2) if not x1 - x2 == 0 else 0
    c = y2 - (m * x2)

    # Using quadratic equation provided in linked solution:
    A = (rad_y ** 2) + ((rad_x ** 2) * (m ** 2))
    B = (2 * (rad_x ** 2) * m * (c - center_y)) - 2 * (rad_y ** 2) * center_x
    C = (rad_y ** 2) * (center_x ** 2) + (rad_x ** 2) * (c - center_y) * 2 - (rad_x ** 2) * (rad_y ** 2)

    # Solving quadratic equation:
    xi1 = (-B - sqrt((B ** 2) - (4 * A * C))) / (2 * A)
    xi2 = (-B + sqrt((B ** 2) - (4 * A * C))) / (2 * A)

    # Solving for y:
    yi1 = m * xi1 + c
    yi2 = m * xi2 + c

    return [[xi1, yi1], [xi2, yi2]]

How do I solve the equation and formula correctly in python without using scipy, numpy etc. (or trigonometry and other solutions described on stackoverflow), in purely mathematical fashion, to retrieve the correct intersection points that lie on the ellipse?

See the explanations in the code how the intersection points are calculated. The code does not use any modules and does all the calculations in pure Python:

def get_intersection_points(line_point_a, line_point_b, center_x, center_y, rad_x, rad_y):
    # Line endpoints
    x1, y1 = line_point_a
    x2, y2 = line_point_b
    # Calculate line slope (m) and intercept (c)
    if x1 != x2:
        m = (y2 - y1) / (x2 - x1)
    else:
        m = float('inf')  # Vertical line case
    c = y1 - m * x1 if m != float('inf') else x1  # If line is vertical, c is x1
    # Ellipse parameters
    h = center_x
    k = center_y
    a = rad_x
    b = rad_y
    if m != float('inf'):
        # Coefficients of the quadratic equation
        A = b**2 + a**2 * m**2
        B = 2 * a**2 * m * (c - k) - 2 * b**2 * h
        C = b**2 * h**2 + a**2 * (c - k)**2 - a**2 * b**2
        # Discriminant
        D = B**2 - 4 * A * C
        if D < 0:
            return []  # No intersection points
        # Calculate the x-coordinates of the intersection points
        sqrt_D = D**0.5 # math.sqrt(D)
        xi1 = (-B - sqrt_D) / (2 * A)
        xi2 = (-B + sqrt_D) / (2 * A)
        # Calculate the y-coordinates of the intersection points
        yi1 = m * xi1 + c
        yi2 = m * xi2 + c
        return [[xi1, yi1], [xi2, yi2]]
    else:
        # Vertical line case (x = c)
        x = c
        A = 1 / a**2
        B = -2 * k / b**2
        C = (h**2 / a**2) + (k**2 / b**2) - 1
        # Quadratic equation for y
        A_y = 1 / b**2
        B_y = -2 * k / b**2
        C_y = (x - h)**2 / a**2 + k**2 / b**2 - 1
        # Discriminant
        D_y = B_y**2 - 4 * A_y * C_y
        
        if D_y < 0:
            return []  # No intersection points
        
        # Calculate the y-coordinates of the intersection points
        sqrt_D_y = D_y**0.5 # math.sqrt(D_y)
        yi1 = (-B_y - sqrt_D_y) / (2 * A_y)
        yi2 = (-B_y + sqrt_D_y) / (2 * A_y)
        
        return [ [x , yi1], [ x , yi2] ]

# Example usage:
line_point_a = (-7,-3)
line_point_b = (7, 3)
center_x = 0
center_y = 0
rad_x = 7
rad_y = 3

print(get_intersection_points(line_point_a, line_point_b, center_x, center_y, rad_x, rad_y))

from turtle import *
import math
def draw_ellipse(center_x, center_y, rad_x, rad_y):
    # Move the turtle to the starting position without drawing
    up()
    goto(center_x + rad_x, center_y)
    down()
    # Number of steps to make the ellipse smooth
    steps = 360
    # Draw the ellipse
    for step in range(steps + 1):
        # Calculate the angle in radians
        angle = math.radians(step)
        
        # Parametric equations for the ellipse
        x = center_x + rad_x * math.cos(angle)
        y = center_y + rad_y * math.sin(angle)
        
        # Move the turtle to the calculated position
        goto(x, y)

from turtle import *
Screen().setworldcoordinates(-10, -10, 10, 10)
Screen().bgcolor("lightgrey")
color("blue")
shape("circle")
pensize(5)
    
draw_ellipse(center_x, center_y, rad_x, rad_y)

up(); goto( center_x, center_y - rad_x ); color("cyan") 
down(); circle(rad_x, steps=36)  
up(); goto( center_x, center_y - rad_y ); 
down(); circle(rad_y, steps=36); color("green")
up(); goto( *line_point_a)
down(); goto(*line_point_b)
up() 
color("red")
for point in get_intersection_points(line_point_a, line_point_b, center_x, center_y, rad_x, rad_y) : 
    goto(*point)
    stamp()
shape("turtle")
color("black")
goto(-8,-8)
done()

To check if the code works correctly the Turtle draws the ellipse, two circles, the line and the found intersection points resting after the done work in the left lower corner of the image: