Intersection of Two Ellipses with a Common Focus, Lua

I am currently working on a project in Lua that involves calculating the intersection points of two ellipses that share a common focus O. The parameters for the ellipses include the coordinates of their shared focus, the coordinates of the second foci (F1 and F2), and the distances from the foci to the vertices of each ellipse (p1 and p2).

(it's a part of Voronoi Diagram, where the sweep line is a circle, the beach line is a fan of ellipses, the project will be written in Lua)

The code for eliipses:

local f0 = {x=400, y=300, r=290}

local function newEllipse (x, y)
    local r = f0.r
    local ellipse = {f0 = f0, f= {x=x, y=y}}

    ellipse.c = length(f0.x, f0.y, x, y)/2
    print ('c', ellipse.c)
    local p = (f0.r - 2*ellipse.c)/2
    print ('p', p)
    ellipse.a = ellipse.c + p
    print ('a', ellipse.a)
    ellipse.b = math.sqrt(ellipse.a^2 - ellipse.c^2)
    
    ellipse.alpha = math.atan2 (y-f0.y, x-f0.x)
    print ('alpha', ellipse.alpha)

    local vertices = {}
    local numPoints = 32
    for i = 1, numPoints do
        local t = (i / numPoints) * (2 * math.pi)
        local x = (f0.x + x) / 2 + ellipse.a * math.cos(t) * math.cos(ellipse.alpha) - ellipse.b * math.sin(t) * math.sin(ellipse.alpha)
        local y = (f0.y + y) / 2 + ellipse.a * math.cos(t) * math.sin(ellipse.alpha) + ellipse.b * math.sin(t) * math.cos(ellipse.alpha)
--      print (x, y)
        table.insert (vertices, x)
        table.insert (vertices, y)
    end
    ellipse.vertices = vertices

    return ellipse
end

local e1 = newEllipse (500, 300)
local e2 = newEllipse (400, 500)

How can I find the points T1 and T2?

Actually, it looks like a solution for PPC Apollonius Problem: Each circle's middle point is the ellipse and the radius to the own focus of it, it touches the circle with tangent!

Update: I've found the solution for not inclined ellipses:

x = f*d
y = f*sqrt(1-d^2)
-- where d is cos theta:
d = (t1 - t2)/(t1*e2 - t2*e1)
t1 = 1-e1^2 -- also t1 = b^2/a^2
t2 = 1-e2^2
-- f is radius for this theta:
f = a*t1/(1-e1*d)

Example:

https://www.desmos.com/calculator/5mhzyskoei

This is almost a purely mathematical problem.

Ellipse centered on a focus in polar coordinates is

ρ = ep / (1 - ecos(θ - Δ))

(ρ, θ) is a point on the ellipse

e is the eccentricity, calculated from c / a

Δ is the rotation of ellipse (the angle from x-axis to the major axis of ellipse)

Now you have 2 ellipses. The goal is to find 2 θ values that satisfy the following equation:

e₁p₁ / (1 - e₁cos(θ - Δ₁)) = e₂p₂ / (1 - e₂cos(θ - Δ₂))

Simplify this equation yields: (You need this formula: cos(a - b) = cos(a)cos(b) + sin(a)sin(b))

(e₁p₁e₂sinΔ₂ - e₂p₂e₁sinΔ₁)sinθ + (e₁p₁e₂cosΔ₂ - e₂p₂e₁cosΔ₁)cosθ = e₁p₁ - e₂p₂

The simplified equation takes the form of Asinθ + Bcosθ = C, so we can use another formula

Asinθ + Bcosθ = sqrt(A² + B²)sin(θ + φ)
φ = atan(b, a)

to calculate the value of θ.

θ = arcsin(C / sqrt(A² + B²)) - φ

---

Other tips you may need:

1. To calculate the angle (Δ) between 2 vectors:

   angle = arccos(dot(v₁, v₂) / sqrt(v₁² * v₂²))
   

   The direction of rotation (from v₁ to v₂) can be determined by: (>0 means clockwise)

   sign = v₁.x * v₂.y - v₁.y * v₂.x
   

2. Since math.asin only returns values between -π/2 to π/2, and since sin(π-θ) = sin(θ), after getting one θ, the other one can be calculated using π-θ

3. Here are the two points I calculated for verification (with the common focus as the origin)

   rho, theta
   P1: 46.845589263185, 0.51139848585219
   P2: 24.79698629529, 3.5149426136911
   
   x, y
   P1: 40.852208863237, 22.926104431619
   P2: -23.088739743658, -9.044369871829