Draw ellipse in matplotlib given the focii

Question

The equation of an ellipse is:

sqrt((x-a1)**2 + (y-b1)**2) + np.sqrt((x-a2)**2 + (y-b2)**2) = c

The focii are (a1, b1) and (a2, b2). c is also known. How do I draw this in python using matplotlib?

Thanks for your help.

Answer 1

You can represent the ellipse parametrically in some variable t. You could look into Wikipedia to see how this can be done, for instance.

In the following code, I've derived the parameters necessary for the parametric form from the parameters that you've supplied.

# Example focii and sum-distance
a1 = 1
b1 = 2
a2 = 5
b2 = 7
c = 9

# Compute ellipse parameters
a = c / 2                                # Semimajor axis
x0 = (a1 + a2) / 2                       # Center x-value
y0 = (b1 + b2) / 2                       # Center y-value
f = np.sqrt((a1 - x0)**2 + (b1 - y0)**2) # Distance from center to focus
b = np.sqrt(a**2 - f**2)                 # Semiminor axis
phi = np.arctan2((b2 - b1), (a2 - a1))   # Angle betw major axis and x-axis

# Parametric plot in t
resolution = 1000
t = np.linspace(0, 2*np.pi, resolution)
x = x0 + a * np.cos(t) * np.cos(phi) - b * np.sin(t) * np.sin(phi)
y = y0 + a * np.cos(t) * np.sin(phi) + b * np.sin(t) * np.cos(phi)

# Plot ellipse
plt.plot(x, y)

# Show focii
plt.plot(a1, b1, 'bo')
plt.plot(a2, b2, 'bo')

plt.axis('equal')
plt.show()

This gives what you need:

Answer 2

you need 2 lists or arrays of X, Y such that the elements satisfy the ellipse equation

the usual ellipse plotting solutions parameterize the ellipse equation by central (or focal) angle to make the X,Y functions single valued for the angle in 0 to 2pi

I showed a solution in Drawing elliptical orbit in Python (using numpy, matplotlib) Y as a function of X with a hack to "feel out" the xrange, then piece together the dual Y solutions for each x

just the minimum mods to that code to put in your equation, would fail for a1 = a2

the symbolic solution takes a minute or so runtime

import numpy as np
import matplotlib.pyplot as plt
from sympy import *  

# sqrt((x-a1)**2 + (y-b1)**2) + np.sqrt((x-a2)**2 + (y-b2)**2) = c
coeffs = [1, 0, -1, 0, 4]
xs = [coeffs[0], coeffs[2]]

def ysolv(coeffs):
    x,y,a1,b1,a2,b2,c = symbols('x y a1 b1 a2 b2 c', real = True)
    ellipse = sqrt((x-a1)**2 + (y-b1)**2) + sqrt((x-a2)**2 + (y-b2)**2) - c
    y_sols = solve(ellipse, y)
    print(*y_sols, sep='\n')

    num_coefs = [(a, f) for a, f in (zip([a1,b1,a2,b2,c], coeffs))]
    y_solsf0 = y_sols[0].subs(num_coefs)
    y_solsf1 = y_sols[1].subs(num_coefs)
    print(y_solsf0, '\n', y_solsf1)

    f0 = lambdify([x], y_solsf0)
    f1 = lambdify([x], y_solsf1)
    return f0, f1

f0, f1 = ysolv(coeffs)

y0 = [f0(x) for x in xs]
y1 = [f1(x) for x in xs]

def feeloutXrange(f, midx, endx):
    fxs = []
    x = midx
    while True:
        try: f(x)
        except:
            break
        fxs.append(x)
        x += (endx - midx)/200
    return fxs

midx = (min(xs) + max(xs))/2    

xpos = feeloutXrange(f0, midx, max(xs))
xnegs = feeloutXrange(f0, midx, min(xs))
xs_ellipse = xnegs[::-1] + xpos[1:]

y0s = [f0(x) for x in xs_ellipse]
y1s = [f1(x) for x in xs_ellipse]

ys_ellipse = y0s + y1s[::-1] + [y0s[0]] # add y start point to end to close drawing

xs_ellipse = xs_ellipse + xs_ellipse[::-1] + [xs_ellipse[0]] # added x start point  

plt.plot(xs_ellipse, ys_ellipse)
plt.show()

(-c*sqrt((a1**2 - 2*a1*a2 + a2**2 + b1**2 - 2*b1*b2 + b2**2 - c**2)*(a1**2 + 2*a1*a2 - 4*a1*x + a2**2 - 4*a2*x + b1**2 - 2*b1*b2 + b2**2 - c**2 + 4*x**2))*(-b1 + b2 + c)*(b1 - b2 + c) + (b1**2 - 2*b1*b2 + b2**2 - c**2)*(-a1**2*b1 + a1**2*b2 + 2*a1*b1*x - 2*a1*b2*x + a2**2*b1 - a2**2*b2 - 2*a2*b1*x + 2*a2*b2*x - b1**3 + b1**2*b2 + b1*b2**2 + b1*c**2 - b2**3 + b2*c**2))/(2*(-b1 + b2 + c)*(b1 - b2 + c)*(b1**2 - 2*b1*b2 + b2**2 - c**2))
(c*sqrt((a1**2 - 2*a1*a2 + a2**2 + b1**2 - 2*b1*b2 + b2**2 - c**2)*(a1**2 + 2*a1*a2 - 4*a1*x + a2**2 - 4*a2*x + b1**2 - 2*b1*b2 + b2**2 - c**2 + 4*x**2))*(-b1 + b2 + c)*(b1 - b2 + c) + (b1**2 - 2*b1*b2 + b2**2 - c**2)*(-a1**2*b1 + a1**2*b2 + 2*a1*b1*x - 2*a1*b2*x + a2**2*b1 - a2**2*b2 - 2*a2*b1*x + 2*a2*b2*x - b1**3 + b1**2*b2 + b1*b2**2 + b1*c**2 - b2**3 + b2*c**2))/(2*(-b1 + b2 + c)*(b1 - b2 + c)*(b1**2 - 2*b1*b2 + b2**2 - c**2))
sqrt(-48*x**2 + 192)/8 
 -sqrt(-48*x**2 + 192)/8

the other answers used the parametric transform approach

I especially like showing sympy solving the equation for you rather than someone hand solving it

the symbolic expression only needs be found once for a particular Ellipse parameterization, then the symbolic expression can simply be hard coded:

"""   
for Ellipse equation:
sqrt((x-a1)**2 + (y-b1)**2) + sqrt((x-a2)**2 + (y-b2)**2) = c   

sympy solution to Ellipse equation, only have to run once to get y_sols
symbolic expression to paste into ysolv below

#def symEllipse():
#    x,y,a1,b1,a2,b2,c = symbols('x y a1 b1 a2 b2 c', real = True)
#    ellipse = sqrt((x-a1)**2 + (y-b1)**2) + sqrt((x-a2)**2 + (y-b2)**2) - c
#    y_sols = solve(ellipse, y)
#    print(*y_sols, sep='\n')

"""

coeffs = [1, 1, -1, -1, 3]
xs = [coeffs[0], coeffs[2]]

def ysolv(coeffs):

    x,y,a1,b1,a2,b2,c = symbols('x y a1 b1 a2 b2 c', real = True)

    y_sols = [
        (-c*sqrt((a1**2 - 2*a1*a2 + a2**2 + b1**2 - 2*b1*b2 + b2**2 - c**2)*
          (a1**2 + 2*a1*a2 - 4*a1*x + a2**2 - 4*a2*x + b1**2 - 2*b1*b2 + b2**2
           - c**2 + 4*x**2))*(-b1 + b2 + c)*(b1 - b2 + c) + (b1**2 - 2*b1*b2 +
           b2**2 - c**2)*(-a1**2*b1 + a1**2*b2 + 2*a1*b1*x - 2*a1*b2*x +
           a2**2*b1 - a2**2*b2 - 2*a2*b1*x + 2*a2*b2*x - b1**3 + b1**2*b2 +
           b1*b2**2 + b1*c**2 - b2**3 + b2*c**2))/(2*(-b1 + b2 + c)*
           (b1 - b2 + c)*(b1**2 - 2*b1*b2 + b2**2 - c**2)),
          (c*sqrt((a1**2 - 2*a1*a2 + a2**2 + b1**2 - 2*b1*b2 + b2**2 - c**2)*
          (a1**2 + 2*a1*a2 - 4*a1*x + a2**2 - 4*a2*x + b1**2 - 2*b1*b2 + b2**2
           - c**2 + 4*x**2))*(-b1 + b2 + c)*(b1 - b2 + c) + (b1**2 - 2*b1*b2 +
           b2**2 - c**2)*(-a1**2*b1 + a1**2*b2 + 2*a1*b1*x - 2*a1*b2*x + 
           a2**2*b1 - a2**2*b2 - 2*a2*b1*x + 2*a2*b2*x - b1**3 + b1**2*b2 + 
           b1*b2**2 + b1*c**2 - b2**3 + b2*c**2))/(2*(-b1 + b2 + c)*
           (b1 - b2 + c)*(b1**2 - 2*b1*b2 + b2**2 - c**2))
          ]
    num_coefs = [(a, f) for a, f in (zip([a1,b1,a2,b2,c], coeffs))]
    y_solsf0 = y_sols[0].subs(num_coefs)
    y_solsf1 = y_sols[1].subs(num_coefs)
    print(y_solsf0, '\n', y_solsf1)

    f0 = lambdify([x], y_solsf0)
    f1 = lambdify([x], y_solsf1)
    return f0, f1