Skip to content

Instantly share code, notes, and snippets.

@PM2Ring
Created November 2, 2023 04:17
Show Gist options
  • Save PM2Ring/a700ff27524c77e945e33040afd8b591 to your computer and use it in GitHub Desktop.
Save PM2Ring/a700ff27524c77e945e33040afd8b591 to your computer and use it in GitHub Desktop.
Plot filled ellipse in 3D
""" 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')
@PM2Ring
Copy link
Author

PM2Ring commented Nov 2, 2023

Live demo on SageMathCell.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment