Plot binary orbits, with either body or the barycentre at the origin

BinaryOrbit.py

""" Kepler binary orbit plot, using eccentric anomaly
    Written by PM 2Ring 2023.10.05
"""

from functools import lru_cache

@lru_cache(int(3))
def kepler(ec, N):
    """ Solve Kepler's equation for N equal timesteps """
    out = []
    E = M = 0
    dM = 2 * pi / N
    for i in range(N):
        M += dM
        E += n(dM / (1 - ec * cos(E)))
        for j in range(9):
            dE = n((M - E + ec * sin(E)) / (1 - ec * cos(E)))
            E += dE
            if abs(dE) < 1e-12:
                break
        out.append(E)
    return out

var('t')
@interact
def orbit(ec=1/5, m1=1/3,
  origin=Selector(["Bary", "Body 1", "Body 2"], selector_type='radio'),
  N=24, size=5, auto_update=False):
    colors = [hue(i/N) for i in srange(N)]

    # Major & minor axes, and focal length
    a = 1
    f = a * ec
    b = sqrt(a^2 - f^2)

    # Total mass is 1
    m2 = 1 - m1

    # Orbit vector function using eccentric anomaly
    xy = vector((a * cos(t) - f, b * sin(t)))

    # Radial vectors between bodies
    # Body vectors obey r1*m1 + r2*m2 = 0 and r1 - r2 = R
    R = [xy(t=E) for E in kepler(ec, N)]

    # Origin
    P = point((0, 0), size=20, color="black", zorder=5, axes=False, figsize=size)

    if origin == "Bary":
        q1, q2 = m2, -m1
        l1, l2 = ":", "--"
    elif origin == "Body 1":
        q1, q2 = -m2, -1
        l1, l2 = "-", "--"
    elif origin == "Body 2":
        q1, q2 = m1, 1
        l1, l2 = "-", ":"

    p1 = [q1 * r for r in R]
    p2 = [q2 * r for r in R]

    P += parametric_plot(xy * q1, (t, 0, 2*pi), color="#aaa", linestyle=l1)
    P += sum(point(r, color=c) for r, c in zip(p1, colors))

    P += parametric_plot(xy * q2, (t, 0, 2*pi), color="#aaa", linestyle=l2)
    P += sum(point(r, color=c) for r, c in zip(p2, colors))

    # Lines connecting the bodies
    P += sum(line([r1, r2], color=c, thickness=0.5) for r1, r2, c in zip(p1, p2, colors))

    P.show()

m1 is the mass of the first body, 0 < m1 < 1, total mass is 1.

Live version on SageMathCell