Earth Moon Distance
From wikiluntti
Introduction
The code
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
# ============================================================
# PHYSICAL PARAMETERS
# ============================================================
G = 6.67430e-11
M_EARTH = 5.972e24
# Real Earth-Moon distance
R = 384_400e3 # metres
# Approximate circular orbital velocity
V = 1022.0 # m/s
DISTANCE_SCALE = 1e9 # metres -> million km
# ============================================================
# VISUAL PARAMETERS
# ============================================================
EARTH_DRAW_RADIUS = 0.064
VECTOR_SCALE0 = 0.15
VECTOR_EVERY = 1
# ============================================================
# ORBITAL PERIOD
# ============================================================
ORBIT_PERIOD = 2 * np.pi * R / V
# ============================================================
# GEOMETRIC CONSTRUCTION
# ============================================================
def calculate_orbit(steps):
# --------------------------------------------------------
# One complete orbit is divided into equal angular steps.
# --------------------------------------------------------
dtheta = 2 * np.pi / steps
tangential_distance = R * np.tan(dtheta)
DT = tangential_distance / V
# Start Moon at x = R
moon = np.array([R, 0.0])
starting_positions = []
tangent_positions = []
moon_positions = []
velocity_vectors = []
gravity_vectors = []
for step in range(steps):
moon_before = moon.copy()
starting_positions.append( moon_before.copy() )
radial = ( moon_before/ np.linalg.norm(moon_before) )
tangent = np.array([-radial[1], radial[0] ])
velocity_vector = tangent * V * DT
velocity_vectors.append( velocity_vector.copy() )
tangent_moon = ( moon_before + velocity_vector )
tangent_positions.append( tangent_moon.copy() )
tangent_distance = np.linalg.norm( tangent_moon )
correction_distance = ( tangent_distance - R )
gravity_direction = ( -tangent_moon / tangent_distance )
gravity_vector = ( gravity_direction * correction_distance )
gravity_vectors.append( gravity_vector.copy() )
moon = (tangent_moon + gravity_vector )
moon_positions.append( moon.copy() )
return (
np.array(starting_positions),
np.array(tangent_positions),
np.array(moon_positions),
np.array(velocity_vectors),
np.array(gravity_vectors),
DT
)
# ============================================================
# INITIAL VALUES
# ============================================================
STEPS0 = 10
(
starting_positions,
tangent_positions,
moon_positions,
velocity_vectors,
gravity_vectors,
DT
) = calculate_orbit(STEPS0)
# ============================================================
# FIGURE
# ============================================================
fig, ax = plt.subplots( figsize=(10, 10) )
plt.subplots_adjust( bottom=0.20)
ax.set_aspect("equal")
ax.set_xlabel("x [million km]")
ax.set_ylabel("y [million km]")
ax.grid( True, alpha=0.3)
ax.set_title( "Geometric construction of one complete Moon orbit")
earth_circle = plt.Circle( (0, 0), EARTH_DRAW_RADIUS, color="royalblue", alpha=0.9, zorder=10)
ax.add_patch( earth_circle)
theta = np.linspace( 0, 2 * np.pi, 500)
orbit_line, = ax.plot(
R * np.cos(theta) / DISTANCE_SCALE,
R * np.sin(theta) / DISTANCE_SCALE,
"--",
color="gray",
linewidth=1,
alpha=0.6,
label="Circular orbit"
)
# ============================================================
# MOON POSITIONS
# ============================================================
moon_plot, = ax.plot(
moon_positions[:, 0] / DISTANCE_SCALE,
moon_positions[:, 1] / DISTANCE_SCALE,
"o",
color="lightgray",
markersize=8,
alpha=0.8,
label="Moon positions"
)
# Highlight starting Moon
moon_initial_plot, = ax.plot(
[starting_positions[0, 0] / DISTANCE_SCALE],
[starting_positions[0, 1] / DISTANCE_SCALE],
"o",
color="dimgray",
markersize=11,
zorder=20
)
# ============================================================
# VECTOR INDICES
# ============================================================
indices = np.arange( 0, STEPS0, VECTOR_EVERY )
velocity_positions = ( starting_positions[indices] / DISTANCE_SCALE)
velocity_u = ( velocity_vectors[indices, 0] / DISTANCE_SCALE * VECTOR_SCALE0 )
velocity_v = ( velocity_vectors[indices, 1] / DISTANCE_SCALE * VECTOR_SCALE0 )
velocity_quiver = ax.quiver(
velocity_positions[:, 0],
velocity_positions[:, 1],
velocity_u,
velocity_v,
color="green",
angles="xy",
scale_units="xy",
scale=1,
width=0.004,
alpha=0.7,
label="Tangential velocity"
)
gravity_positions = ( tangent_positions[indices] / DISTANCE_SCALE)
gravity_u = ( gravity_vectors[indices, 0] / DISTANCE_SCALE * VECTOR_SCALE0)
gravity_v = ( gravity_vectors[indices, 1] / DISTANCE_SCALE * VECTOR_SCALE0)
gravity_quiver = ax.quiver(
gravity_positions[:, 0],
gravity_positions[:, 1],
gravity_u,
gravity_v,
color="red",
angles="xy",
scale_units="xy",
scale=1,
width=0.004,
alpha=0.7,
label="Gravity / correction"
)
info = ax.text(
0.02,
0.98,
"",
transform=ax.transAxes,
verticalalignment="top",
fontsize=10,
bbox=dict(
boxstyle="round",
facecolor="white",
alpha=0.85
)
)
# Steps slider
ax_steps = plt.axes( [0.15, 0.10, 0.70, 0.03])
slider_steps = Slider(
ax_steps,
"Steps / orbit",
5,
200,
valinit=STEPS0,
valstep=1
)
# Vector scale slider
ax_scale = plt.axes( [0.15, 0.05, 0.70, 0.03])
slider_scale = Slider(
ax_scale,
"Vector scale",
1,
1.0,
valinit=VECTOR_SCALE0
)
# ============================================================
# UPDATE FUNCTION
# ============================================================
def update(_):
global velocity_quiver, gravity_quiver
steps = int( slider_steps.val)
vector_scale = ( slider_scale.val)
(
starting_positions,
tangent_positions,
moon_positions,
velocity_vectors,
gravity_vectors,
DT
) = calculate_orbit(steps)
moon_plot.set_data( moon_positions[:, 0] / DISTANCE_SCALE, moon_positions[:, 1] / DISTANCE_SCALE )
# Starting Moon position
moon_initial_plot.set_data([starting_positions[0, 0] / DISTANCE_SCALE],[starting_positions[0, 1] / DISTANCE_SCALE] )
indices = np.arange( 0, steps, VECTOR_EVERY)
velocity_quiver.remove()
gravity_quiver.remove()
velocity_positions = ( starting_positions[indices] / DISTANCE_SCALE )
velocity_u = ( velocity_vectors[indices, 0] / DISTANCE_SCALE * vector_scale )
velocity_v = ( velocity_vectors[indices, 1] / DISTANCE_SCALE * vector_scale )
velocity_quiver = ax.quiver(
velocity_positions[:, 0],
velocity_positions[:, 1],
velocity_u,
velocity_v,
color="green",
angles="xy",
scale_units="xy",
scale=1,
width=0.004,
alpha=0.7
)
gravity_positions = ( tangent_positions[indices] / DISTANCE_SCALE )
gravity_u = ( gravity_vectors[indices, 0] / DISTANCE_SCALE * vector_scale)
gravity_v = ( gravity_vectors[indices, 1] / DISTANCE_SCALE * vector_scale)
gravity_quiver = ax.quiver(
gravity_positions[:, 0],
gravity_positions[:, 1],
gravity_u,
gravity_v,
color="red",
angles="xy",
scale_units="xy",
scale=1,
width=0.004,
alpha=0.7
)
green_length = np.linalg.norm(velocity_vectors[0])
red_length = np.linalg.norm(gravity_vectors[0])
green_fraction = (green_length / R * 100)
red_fraction = (red_length / R * 100)
# ========================================================
# INFORMATION
# ========================================================
circular_velocity = np.sqrt(G * M_EARTH / R)
info.set_text(
f"Earth-Moon distance = {R / 1000:,.0f} km\n"
f"Moon velocity = {V:.1f} m/s\n"
f"Circular velocity = {circular_velocity:.1f} m/s\n"
f"Orbit period = {ORBIT_PERIOD / 86400:.2f} days\n"
f"Steps / orbit = {steps}\n"
f"DT per step = {DT / 86400:.2f} days\n"
f"\n"
f"Green vector = {green_length / 1000:,.0f} km"
f" ({green_fraction:.1f}% of R)\n"
f"Red vector = {red_length / 1000:,.0f} km"
f" ({red_fraction:.1f}% of R)"
)
# ========================================================
# REDRAW
# ========================================================
fig.canvas.draw_idle()
slider_steps.on_changed(update)
slider_scale.on_changed(update)
update(None)
ax.legend(loc="lower right")
plt.show()
