Plot filled ellipse in 3D

ellipse3d.py

""" Plot filled ellipse in 3D.
    With major axis and line of nodes
    Written by PM 2Ring 2023.9.24
"""

var('th')

deg = n(pi / 180)

def ellipse(a, ec):
    p = a * (1 - ec^2)
    r = p / (1 + ec * cos(th))
    return vector([r * cos(th), r * sin(th)])

def filled_plot(func, dom, color="#cc0"):
    color = Color(color)
    return parametric_plot(func, dom, fill=True,
      color=color.darker(0.5), fillcolor=color,
      alpha=0.9, fillalpha=0.6, thickness=2
    )

# XY plane & axes
def add_axes(P):
    bbox = P.bounding_box()
    # XY plane
    Q = plot3d(0, *[*zip(*bbox)][:2], color="#ccc", opacity=0.6)

    dlo, dhi = [matrix.diagonal(t) for t in bbox]
    colors = "red", "green", "blue"
    # Axes
    Q += sum(line3d([lo, hi], color=c)
      for lo, hi, c in zip(dlo, dhi, colors))
    Q += sum(point(p, size=5, color=c)
      for p, c in zip(dhi, colors))
    return Q

# Orbital elements. Angles in radians
# a: Semi-major axis
# ec: Eccentricity
# Om: Long of Ascending Node
# In: Inclination
# W: Arg of perifocus

def plot_orbit(a, ec, Om, In, W, color="#cc0"):
    ell = ellipse(a, ec)
    P = filled_plot(ell, (th, 0, 2*pi), color=color)
    # Periapsis, apoapsis
    rp, ra = a * (1 - ec), a * (1 + ec)
    # Major axis
    P += line([(rp, 0), (-ra, 0)], color="#500")
    # Line of nodes
    P += line([ell(th=-W), ell(th=pi-W)], color="#033")
    return P.plot3d().rotateZ(-W).rotateX(-In).rotateZ(-Om)

P = Graphics()

a = 10
ec = 0
In = 30 * deg
W = 0

mi = 4
for i in range(mi):
    Om = i * 30 * deg
    P += plot_orbit(a, ec, Om, In, W, color=hue(i/mi))

P += add_axes(P)
P.show(frame=False, theme='dark', projection='orthographic')

Live demo on SageMathCell.