Skip to content

Instantly share code, notes, and snippets.

@ad-1
Created June 24, 2024 18:47
Show Gist options
  • Select an option

  • Save ad-1/3230201f98907e9b3f04e9ad79992023 to your computer and use it in GitHub Desktop.

Select an option

Save ad-1/3230201f98907e9b3f04e9ad79992023 to your computer and use it in GitHub Desktop.
Elliptical Anomalies
df = pd.DataFrame({
'Position (m)': r,
'Position (X)': rx,
'Position (Y)': ry,
'True Anomaly (rad)': theta,
'True Anomaly (deg)': rad2deg(theta),
'Eccentric Anomaly (rad)': E,
'Eccentric Anomaly (deg)': rad2deg(E),
'Eccentric Anomaly (X)': eccentric_anomaly_x,
'Eccentric Anomaly (Y)': eccentric_anomaly_y,
'Mean Anomaly (rad)': Me,
'Mean Anomaly (deg)': rad2deg(Me),
'Mean Anomaly (X)': mean_anomaly_x,
'Mean Anomaly (Y)': mean_anomaly_y,
'Time Since Periapsis (s)': tp,
'Time Since Periapsis (hr)': tp / (60 * 60),
})
# Create a Matplotlib figure and axis
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_facecolor(colors["background-secondary"])
fig.patch.set_facecolor(colors["background-primary"])
# Plot the orbit and auxiliary circle
ax.plot(df['Position (X)'], df['Position (Y)'], label='Ellipse', color=colors['primary-color'])
true_anomaly_plot_line, = ax.plot([], [], 'p-', label='True Anomaly', color=colors['warning-color'])
true_anomaly_text = ax.annotate("", (c, Re / 2), fontsize=10, ha='left', color=colors["warning-color"])
ax.plot(auxiliary_circle_x, auxiliary_circle_y, label='Auxiliary Circle', color=colors["accent-light"])
eccentric_anomaly_plot_line, = ax.plot([], [], 'p-', label='Eccentric Anomaly', color=colors["accent-light"])
eccentric_anomaly_text = ax.annotate("", (0, Re / 2), fontsize=10, ha='left', color=colors["accent-light"])
eccentric_vertical_plot_line, = ax.plot([], [], 'p-', label='Vertical', color=colors["accent-dark"])
mean_anomaly_plot_line, = ax.plot([], [], 'p-', label='Mean Anomaly', color=colors["info-color"])
mean_anomaly_text = ax.annotate("", (0, -Re), fontsize=10, ha='left', color=colors["info-color"])
time_since_periapsis_text = ax.annotate("", (-a, -a), fontsize=10, ha='left', color=colors["text-primary"])
# Customize the appearance of the text and labels
ax.set_title('Elliptical Orbit Anomalies', color=colors["text-primary"])
ax.set_xlabel('X-axis', color=colors["text-primary"])
ax.set_ylabel('Y-axis', color=colors["text-primary"])
# plot Foci
ax.text(-c, -Re, 'F1', fontsize=12, ha='right', color=colors['text-primary'])
ax.text(c, -Re, 'F2', fontsize=12, ha='left', color=colors['text-primary'])
ax.scatter([-c, c], [0, 0], color=colors['text-primary'], marker='o', label='Foci')
# Function to update the line and point's position during animation
def plot_anomalies_in_time(
frame: int,
) -> Tuple[Line2D, Line2D, Line2D, Annotation, Line2D, Annotation, Annotation, Annotation]:
# Update the position of the true anomaly
true_anomaly_plot_line.set_data(
[c, df['Position (X)'][frame]],
[0, df['Position (Y)'][frame]]
)
# Extract the position of the eccentric anomaly line
eccentric_anomaly_plot_line.set_data(
[0, df['Eccentric Anomaly (X)'][frame]],
[0, df['Eccentric Anomaly (Y)'][frame]]
)
eccentric_vertical_plot_line.set_data(
[eccentric_anomaly_x[frame], eccentric_anomaly_x[frame]],
[df['Position (Y)'][frame], eccentric_anomaly_y[frame]]
)
mean_anomaly_plot_line.set_data(
[0, df['Mean Anomaly (X)'][frame]],
[0, df['Mean Anomaly (Y)'][frame]]
)
eccentric_anomaly_text.set_text(f'E = {df["Eccentric Anomaly (deg)"][frame]:.1f}°')
mean_anomaly_text.set_text(f'Me = {df["Mean Anomaly (deg)"][frame]:.1f}°')
true_anomaly_text.set_text(rf'$\theta$ = {df["True Anomaly (deg)"][frame]:.1f}°')
time_since_periapsis_text.set_text(f'tp = {df["Time Since Periapsis (hr)"][frame]:.2f} hr')
return (
true_anomaly_plot_line,
eccentric_anomaly_plot_line,
eccentric_vertical_plot_line,
eccentric_anomaly_text,
mean_anomaly_plot_line,
mean_anomaly_text,
true_anomaly_text,
time_since_periapsis_text,
)
# Time steps
n_frames = len(theta)
frames = np.arange(0, n_frames)
# Create the animation
ani = FuncAnimation(fig, plot_anomalies_in_time, frames=frames, blit=True)
ani.save('elliptical-anomalies.gif', writer='pillow')
ax.set_aspect('equal')
ax.grid(color=colors["background-card"])
ax.legend(loc='upper right', bbox_to_anchor=(1.3, 1.0))
plt.show()
def calc_angular_momentum(standard_gravitational_parameter, position, eccentricity, true_anomaly):
"""
Calculate the angular momentum of an orbiting body.
Parameters:
standard_gravitational_parameter (float): The standard gravitational parameter (GM) of the primary body.
position (float): The current position (radius) of the orbiting body.
eccentricity (float): The eccentricity of the orbit.
true_anomaly (float): The true anomaly (angle) at the given position.
Returns:
float: The angular momentum of the orbiting body.
"""
return np.sqrt(standard_gravitational_parameter * position * (1 + eccentricity * np.cos(true_anomaly)))
# Orbit eccentricity
e = calc_eccentricity(
radius_apoapsis=ra,
radius_periapsis=rp
)
# Orbit angular momentum
h = calc_angular_momentum(
standard_gravitational_parameter=mu,
position=rp,
eccentricity=e,
true_anomaly=true_anomaly_at_periapsis
)
# Orbit period
T = calc_period_from_momentum(
standard_gravitational_parameter=mu,
angular_momentum=h,
eccentricity=e
)
# True Anomaly
num_points = 180 # Number of angle points
theta = deg2rad(np.linspace(0, 360, num_points))
# Eccentric Anomaly
E = calc_eccentric_anomaly(eccentricity=e, true_anomaly=theta)
# Mean anomaly
Me = calc_mean_anomaly(eccentric_anomaly=E, eccentricity=e)
# Time since periapsis
tp = calc_time_since_periapsis(mean_anomaly=Me, period=T)
# Elliptical geometry
a = calc_semi_major_axis(periapsis_radius=rp, apoapsis_radius=ra)
c = calc_foci_position(semi_major_axis=a, eccentricity=e)
b = calc_semi_minor_axis(semi_major_axis=a, eccentricity=e)
# Orbit position
r = orbit_position(
angular_momentum=h,
eccentricity=e,
standard_gravitational_parameter=mu,
true_anomaly=theta
)
def calc_eccentric_anomaly(eccentricity, true_anomaly):
"""
Calculate the eccentric anomaly from the true anomaly and eccentricity.
Parameters:
eccentricity (float): The eccentricity of the orbit.
true_anomaly (float): The true anomaly.
Returns:
np.array: The eccentric anomaly.
"""
_E = 2 * (np.arctan(np.sqrt((1 - eccentricity)/(1 + eccentricity)) * np.tan(true_anomaly / 2)))
return np.array([value if value >= 0 else 2 * np.pi - abs(value) for value in _E])
def calc_eccentricity(radius_apoapsis, radius_periapsis):
"""
Calculate the eccentricity of an orbit.
Parameters:
radius_apoapsis (float): The radius of the apoapsis (farthest point in orbit).
radius_periapsis (float): The radius of the periapsis (closest point in orbit).
Returns:
float: The eccentricity of the orbit.
"""
return (radius_apoapsis - radius_periapsis) / (radius_apoapsis + radius_periapsis)
def calc_foci_position(semi_major_axis, eccentricity):
"""
Calculate the position of the foci of the ellipse.
Parameters:
semi_major_axis (float): The semi-major axis of the orbit.
eccentricity (float): The eccentricity of the orbit.
Returns:
float: The distance from the center of the ellipse to each focus.
"""
return semi_major_axis * eccentricity
# Calculate the x and y coordinates of the orbit in Cartesian Space
rx = (a * e) + r * np.cos(theta)
ry = np.multiply(r, np.sin(theta))
# Calculate the x and y coordinates of the auxiliary circle
auxiliary_circle_x = a * np.cos(theta)
auxiliary_circle_y = a * np.sin(theta)
# Calculate the x and y coordinates of the eccentric anomaly
eccentric_anomaly_x = a * np.cos(E)
eccentric_anomaly_y = a * np.sin(E)
# Calculate the x and y coordinates of the mean anomaly
mean_anomaly_x = a * np.cos(Me)
mean_anomaly_y = a * np.sin(Me)
def calc_mean_anomaly(eccentric_anomaly, eccentricity):
"""
Calculate the mean anomaly from the eccentric anomaly and eccentricity.
Parameters:
eccentric_anomaly (float): The eccentric anomaly.
eccentricity (float): The eccentricity of the orbit.
Returns:
float: The mean anomaly.
"""
return eccentric_anomaly - (eccentricity * np.sin(eccentric_anomaly))
def calc_mean_rate(period):
"""
Calculate the mean motion (rate) of the orbiting body.
Parameters:
period (float): The orbital period.
Returns:
float: The mean motion (rate) of the orbiting body.
"""
return 2 * np.pi / period
def calc_period_from_axis(standard_gravitational_parameter, semi_major_axis):
"""
Calculate the orbital period from the semi-major axis.
Parameters:
standard_gravitational_parameter (float): The standard gravitational parameter (GM) of the primary body.
semi_major_axis (float): The semi-major axis of the orbit.
Returns:
float: The orbital period.
"""
return ((2 * np.pi) / (np.sqrt(standard_gravitational_parameter))) * semi_major_axis ** (3/2)
def calc_period_from_momentum(standard_gravitational_parameter, angular_momentum, eccentricity):
"""
Calculate the orbital period from angular momentum and eccentricity.
Parameters:
standard_gravitational_parameter (float): The standard gravitational parameter (GM) of the primary body.
angular_momentum (float): The angular momentum of the orbiting body.
eccentricity (float): The eccentricity of the orbit.
Returns:
float: The orbital period.
"""
return ((2 * np.pi) / (standard_gravitational_parameter ** 2)) * \
(angular_momentum / np.sqrt(1 - eccentricity ** 2)) ** 3
def calc_semi_major_axis(periapsis_radius, apoapsis_radius):
"""
Calculate the semi-major axis of the orbit.
Parameters:
periapsis_radius (float): The radius of the periapsis (closest point in orbit).
apoapsis_radius (float): The radius of the apoapsis (farthest point in orbit).
Returns:
float: The semi-major axis.
"""
return (periapsis_radius + apoapsis_radius) / 2.0
def calc_semi_minor_axis(semi_major_axis, eccentricity):
"""
Calculate the semi-minor axis of the orbit.
Parameters:
semi_major_axis (float): The semi-major axis of the orbit.
eccentricity (float): The eccentricity of the orbit.
Returns:
float: The semi-minor axis.
"""
return semi_major_axis * np.sqrt(1 - (eccentricity ** 2))
def calc_time_since_periapsis(mean_anomaly, period):
"""
Calculate the time since periapsis from the mean anomaly and orbital period.
Parameters:
mean_anomaly (float): The mean anomaly.
period (float): The orbital period.
Returns:
float: The time since periapsis.
"""
return (mean_anomaly / (2 * np.pi)) * period
# Orbit configuration
Re = 6378
mu = 398600
rp = Re + 9600
ra = Re + 100000
true_anomaly_at_periapsis = 0
def orbit_position(angular_momentum, standard_gravitational_parameter, eccentricity, true_anomaly):
"""
Calculate the position of the orbiting body in its orbit.
Parameters:
angular_momentum (float): The angular momentum of the orbiting body.
standard_gravitational_parameter (float): The standard gravitational parameter (GM) of the primary body.
eccentricity (float): The eccentricity of the orbit.
true_anomaly (float): The true anomaly (angle) at the given position.
Returns:
float: The position (radius) of the orbiting body.
"""
return ((angular_momentum ** 2) / standard_gravitational_parameter) / (1 + (eccentricity * np.cos(true_anomaly)))
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.lines import Line2D
from matplotlib.text import Annotation
from typing import Tuple
colors = {
"primary-color": "#1976D2", # Material Blue 500
"background-primary": "#1E1E1E", # Dark Background
"text-primary": "#FFFFFF", # White text on dark background
"background-secondary": "#2C2C2C", # Slightly lighter background
"background-card": "#3A3A3A", # Grid background
"accent-dark": "#D84315", # Material Deep Orange 600
"accent-light": "#FF8A65", # Material Deep Orange 300
"warning-color": "#FBC02D", # Material Yellow 700 for warnings
"info-color": "#0288D1", # Material Light Blue 600 for informational messages
}
def deg2rad(deg):
return deg * np.pi / 180
def rad2deg(rad):
return rad * 180 / np.pi
def calc_eccentricity(radius_apoapsis, radius_periapsis):
return (radius_apoapsis - radius_periapsis) / (radius_apoapsis + radius_periapsis)
def calc_angular_momentum(standard_gravitational_parameter, position, eccentricity, true_anomaly):
return np.sqrt(standard_gravitational_parameter * position * (1 + eccentricity * np.cos(true_anomaly)))
def calc_period_from_momentum(standard_gravitational_parameter, angular_momentum, eccentricity):
return ((2 * np.pi) / (standard_gravitational_parameter ** 2)) * \
(angular_momentum / np.sqrt(1 - eccentricity ** 2)) ** 3
def calc_period_from_axis(standard_gravitational_parameter, semi_major_axis):
return ((2 * np.pi) / (np.sqrt(standard_gravitational_parameter))) * semi_major_axis ** (3/2)
def calc_mean_anomaly(eccentric_anomaly, eccentricity):
return eccentric_anomaly - (eccentricity * np.sin(eccentric_anomaly))
def calc_time_since_periapsis(mean_anomaly, period):
return (mean_anomaly / (2 * np.pi)) * period
def calc_eccentric_anomaly(eccentricity, true_anomaly):
_E = 2 * (np.arctan(np.sqrt((1 - eccentricity)/(1 + eccentricity)) * np.tan(true_anomaly / 2)))
return np.array([value if value >= 0 else 2 * np.pi - abs(value) for value in _E])
def calc_semi_major_axis(periapsis_radius, apoapsis_radius):
return (periapsis_radius + apoapsis_radius) / 2.0
def calc_semi_minor_axis(semi_major_axis, eccentricity):
return semi_major_axis * np.sqrt(1 - (eccentricity ** 2))
def calc_foci_position(semi_major_axis, eccentricity):
return semi_major_axis * eccentricity
def orbit_position(angular_momentum, standard_gravitational_parameter, eccentricity, true_anomaly):
return ((angular_momentum ** 2) / standard_gravitational_parameter) / (1 + (eccentricity * np.cos(true_anomaly)))
def calc_mean_rate(period):
return 2 * np.pi / period
# Orbit configuration
Re = 6378
mu = 398600
rp = Re + 9600
ra = Re + 100000
true_anomaly_at_periapsis = 0
# Orbit eccentricity
e = calc_eccentricity(
radius_apoapsis=ra,
radius_periapsis=rp
)
# Orbit angular momentum
h = calc_angular_momentum(
standard_gravitational_parameter=mu,
position=rp,
eccentricity=e,
true_anomaly=true_anomaly_at_periapsis
)
# Orbit period
T = calc_period_from_momentum(
standard_gravitational_parameter=mu,
angular_momentum=h,
eccentricity=e
)
# True Anomaly
num_points = 180 # Number of angle points
theta = deg2rad(np.linspace(0, 360, num_points))
# Eccentric Anomaly
E = calc_eccentric_anomaly(eccentricity=e, true_anomaly=theta)
# Mean anomaly
Me = calc_mean_anomaly(eccentric_anomaly=E, eccentricity=e)
# Time since periapsis
tp = calc_time_since_periapsis(mean_anomaly=Me, period=T)
# Elliptical geometry
a = calc_semi_major_axis(periapsis_radius=rp, apoapsis_radius=ra)
c = calc_foci_position(semi_major_axis=a, eccentricity=e)
b = calc_semi_minor_axis(semi_major_axis=a, eccentricity=e)
# Orbit position
r = orbit_position(
angular_momentum=h,
eccentricity=e,
standard_gravitational_parameter=mu,
true_anomaly=theta
)
# Calculate the x and y coordinates of the orbit in Cartesian Space
rx = (a * e) + r * np.cos(theta)
ry = np.multiply(r, np.sin(theta))
# Calculate the x and y coordinates of the auxiliary circle
auxiliary_circle_x = a * np.cos(theta)
auxiliary_circle_y = a * np.sin(theta)
# Calculate the x and y coordinates of the eccentric anomaly
eccentric_anomaly_x = a * np.cos(E)
eccentric_anomaly_y = a * np.sin(E)
# Calculate the x and y coordinates of the mean anomaly
mean_anomaly_x = a * np.cos(Me)
mean_anomaly_y = a * np.sin(Me)
df = pd.DataFrame({
'Position (m)': r,
'Position (X)': rx,
'Position (Y)': ry,
'True Anomaly (rad)': theta,
'True Anomaly (deg)': rad2deg(theta),
'Eccentric Anomaly (rad)': E,
'Eccentric Anomaly (deg)': rad2deg(E),
'Eccentric Anomaly (X)': eccentric_anomaly_x,
'Eccentric Anomaly (Y)': eccentric_anomaly_y,
'Mean Anomaly (rad)': Me,
'Mean Anomaly (deg)': rad2deg(Me),
'Mean Anomaly (X)': mean_anomaly_x,
'Mean Anomaly (Y)': mean_anomaly_y,
'Time Since Periapsis (s)': tp,
'Time Since Periapsis (hr)': tp / (60 * 60),
})
# Create a Matplotlib figure and axis
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_facecolor(colors["background-secondary"])
fig.patch.set_facecolor(colors["background-primary"])
# Plot the orbit and auxiliary circle
ax.plot(df['Position (X)'], df['Position (Y)'], label='Ellipse', color=colors['primary-color'])
true_anomaly_plot_line, = ax.plot([], [], 'p-', label='True Anomaly', color=colors['warning-color'])
true_anomaly_text = ax.annotate("", (c, Re / 2), fontsize=10, ha='left', color=colors["warning-color"])
ax.plot(auxiliary_circle_x, auxiliary_circle_y, label='Auxiliary Circle', color=colors["accent-light"])
eccentric_anomaly_plot_line, = ax.plot([], [], 'p-', label='Eccentric Anomaly', color=colors["accent-light"])
eccentric_anomaly_text = ax.annotate("", (0, Re / 2), fontsize=10, ha='left', color=colors["accent-light"])
eccentric_vertical_plot_line, = ax.plot([], [], 'p-', label='Vertical', color=colors["accent-dark"])
mean_anomaly_plot_line, = ax.plot([], [], 'p-', label='Mean Anomaly', color=colors["info-color"])
mean_anomaly_text = ax.annotate("", (0, -Re), fontsize=10, ha='left', color=colors["info-color"])
time_since_periapsis_text = ax.annotate("", (-a, -a), fontsize=10, ha='left', color=colors["text-primary"])
# Customize the appearance of the text and labels
ax.set_title('Elliptical Orbit Anomalies', color=colors["text-primary"])
ax.set_xlabel('X-axis', color=colors["text-primary"])
ax.set_ylabel('Y-axis', color=colors["text-primary"])
# plot Foci
ax.text(-c, -Re, 'F1', fontsize=12, ha='right', color=colors['text-primary'])
ax.text(c, -Re, 'F2', fontsize=12, ha='left', color=colors['text-primary'])
ax.scatter([-c, c], [0, 0], color=colors['text-primary'], marker='o', label='Foci')
# Function to update the line and point's position during animation
def plot_anomalies_in_time(
frame: int,
) -> Tuple[Line2D, Line2D, Line2D, Annotation, Line2D, Annotation, Annotation, Annotation]:
# Update the position of the true anomaly
true_anomaly_plot_line.set_data(
[c, df['Position (X)'][frame]],
[0, df['Position (Y)'][frame]]
)
# Extract the position of the eccentric anomaly line
eccentric_anomaly_plot_line.set_data(
[0, df['Eccentric Anomaly (X)'][frame]],
[0, df['Eccentric Anomaly (Y)'][frame]]
)
eccentric_vertical_plot_line.set_data(
[eccentric_anomaly_x[frame], eccentric_anomaly_x[frame]],
[df['Position (Y)'][frame], eccentric_anomaly_y[frame]]
)
mean_anomaly_plot_line.set_data(
[0, df['Mean Anomaly (X)'][frame]],
[0, df['Mean Anomaly (Y)'][frame]]
)
eccentric_anomaly_text.set_text(f'E = {df["Eccentric Anomaly (deg)"][frame]:.1f}°')
mean_anomaly_text.set_text(f'Me = {df["Mean Anomaly (deg)"][frame]:.1f}°')
true_anomaly_text.set_text(rf'$\theta$ = {df["True Anomaly (deg)"][frame]:.1f}°')
time_since_periapsis_text.set_text(f'tp = {df["Time Since Periapsis (hr)"][frame]:.2f} hr')
return (
true_anomaly_plot_line,
eccentric_anomaly_plot_line,
eccentric_vertical_plot_line,
eccentric_anomaly_text,
mean_anomaly_plot_line,
mean_anomaly_text,
true_anomaly_text,
time_since_periapsis_text,
)
# Time steps
n_frames = len(theta)
frames = np.arange(0, n_frames)
# Create the animation
ani = FuncAnimation(fig, plot_anomalies_in_time, frames=frames, blit=True)
ani.save('elliptical-anomalies.gif', writer='pillow')
ax.set_aspect('equal')
ax.grid(color=colors["background-card"])
ax.legend(loc='upper right', bbox_to_anchor=(1.3, 1.0))
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment