Points of Intersection of an Ellipse and a line after rotating ellipse by angle theta

I want to find points of intersection of an Ellipse and a line after rotating Ellipse by angle theta.

I have written python code to find points of intersection of an Ellipse and a line, but I am not able to figure out how to find intersection point after rotating the Ellipse by theta angle.

def intersactionPoints(a,b,h,k,x1,y1,x2,y2):
    #xi1, yi1, xi2, yi2 <- intersection points
    xi1, yi1, xi2, yi2, aa, bb, cc, m = 0, 0, 0, 0, 0, 0, 0, 0
    if x1 != x2:
        m = (y2 - y1)/(x2 - x1)
        c = y1 - m * x1
        aa = b * b + a * a * m * m
        bb = 2 * a * a * c * m - 2 * a * a * k * m - 2 * h * b * b
        cc = b * b * h * h + a * a * c * c - 2 * a * a * k * c + a * a * k * k - a * a * b * b
    else:
        # vertical line case
        aa = a * a
        bb = -2.0 * k * a * a
        cc = -a * a * b * b + b * b * (x1 - h) * (x1 - h)
    d = bb * bb - 4 * aa * cc
    # intersection points : (xi1,yi1) and (xi2,yi2)
    if d > 0:
        if (x1 != x2):
            xi1 = (-bb + (d**0.5)) / (2 * aa)
            xi2 = (-bb - (d**0.5)) / (2 * aa)
            yi1 = y1 + m * (xi1 - x1)
            yi2 = y1 + m * (xi2 - x1)
        else:
            yi1 = (-bb + (d**0.5)) / (2 * aa)
            yi2 = (-bb - (d**0.5)) / (2 * aa)
            xi1 = x1
            xi2 = x1
    return xi1, yi1, xi2, yi2


if __name__ == "__main__":
    a =  #major axis
    b =  #minor axis
    h =  #center x of ellipse
    k =  #center y of ellipse
    x1 =  #line coordinate x1
    y1 =  #line coordinate y1
    x2 = #line coordinate x2
    y2 = #line coordinate y2
    xi1, yi1, xi2, yi2 = intersactionPoints(a,b,h,k,x1,y1,x2,y2)

If you have tested and ready-to-use solution for axis-aligned ellipse, it is much simpler to transform points defining line into ellipse system, find intersection points, then make reverse transformation.

For origin-centered ellipse and rotation angle theta

 x1' = x1 * Cos(theta) + y1 * Sin(theta)
 y1' = - x1 * Sin(theta) + y1 * Cos(theta)

For ellipse centered at cx, cy:

 x1' = (x1 - cx) * Cos(theta) + (y1 - cy) * Sin(theta)
 y1' = - (x1 - cx) * Sin(theta) + (y1 - cy) * Cos(theta)