7 changed files with 337 additions and 70 deletions
@ -0,0 +1,337 @@
|
||||
#!/usr/bin/env python3 |
||||
""" |
||||
Precalculate moon orbit TOML config from planetary_data.md values. |
||||
|
||||
Converts mean anomaly (M) at J2000 to true anomaly (nu) via Kepler's equation, |
||||
then outputs a complete test TOML config with correct planetary masses, |
||||
eccentricities, inclinations, and orbital elements. |
||||
|
||||
Usage: |
||||
python3 scripts/precalc_moon_orbits.py |
||||
|
||||
Outputs: |
||||
- TOML config to stdout (redirect to tests/test_moon_orbits.toml) |
||||
- Console summary of computed values |
||||
""" |
||||
|
||||
import sys, math |
||||
|
||||
sys.path.insert(0, "scripts") |
||||
from sim_engine import ( |
||||
solve_kepler_elliptical, |
||||
orbital_to_cartesian, |
||||
vmag, |
||||
G, |
||||
normalize_angle, |
||||
OrbitalElements, |
||||
) |
||||
|
||||
# Planetary data from docs/planetary_data.md |
||||
|
||||
SUN_MASS = 1.989e30 |
||||
SUN_RADIUS = 6.96e8 |
||||
|
||||
PLANETS = [ |
||||
{ |
||||
"name": "Venus", |
||||
"mass": 4.87e24, |
||||
"radius": 6.052e6, |
||||
"parent": 0, |
||||
"a_au": 0.723, |
||||
"e": 0.007, |
||||
"inc_deg": 3.39, |
||||
"Omega_deg": 76.68, |
||||
"omega_deg": 54.92, |
||||
"M_deg": 50.38, |
||||
}, |
||||
{ |
||||
"name": "Earth", |
||||
"mass": 5.97e24, |
||||
"radius": 6.378e6, |
||||
"parent": 0, |
||||
"a_au": 1.000, |
||||
"e": 0.017, |
||||
"inc_deg": 0.00, |
||||
"Omega_deg": 0.00, |
||||
"omega_deg": 102.94, |
||||
"M_deg": -2.47, |
||||
}, |
||||
{ |
||||
"name": "Mars", |
||||
"mass": 6.42e23, |
||||
"radius": 3.396e6, |
||||
"parent": 0, |
||||
"a_au": 1.524, |
||||
"e": 0.093, |
||||
"inc_deg": 1.85, |
||||
"Omega_deg": 49.56, |
||||
"omega_deg": 286.50, |
||||
"M_deg": 19.39, |
||||
}, |
||||
{ |
||||
"name": "Jupiter", |
||||
"mass": 1.898e27, |
||||
"radius": 71.492e6, |
||||
"parent": 0, |
||||
"a_au": 5.203, |
||||
"e": 0.049, |
||||
"inc_deg": 1.31, |
||||
"Omega_deg": 100.47, |
||||
"omega_deg": 274.25, |
||||
"M_deg": 19.67, |
||||
}, |
||||
{ |
||||
"name": "Saturn", |
||||
"mass": 5.683e26, |
||||
"radius": 60.268e6, |
||||
"parent": 0, |
||||
"a_au": 9.537, |
||||
"e": 0.057, |
||||
"inc_deg": 2.49, |
||||
"Omega_deg": 113.66, |
||||
"omega_deg": 338.94, |
||||
"M_deg": -42.64, |
||||
}, |
||||
{ |
||||
"name": "Uranus", |
||||
"mass": 8.68e25, |
||||
"radius": 25.559e6, |
||||
"parent": 0, |
||||
"a_au": 19.19, |
||||
"e": 0.046, |
||||
"inc_deg": 0.77, |
||||
"Omega_deg": 74.02, |
||||
"omega_deg": 96.94, |
||||
"M_deg": 142.28, |
||||
}, |
||||
{ |
||||
"name": "Neptune", |
||||
"mass": 1.02e26, |
||||
"radius": 24.764e6, |
||||
"parent": 0, |
||||
"a_au": 30.07, |
||||
"e": 0.010, |
||||
"inc_deg": 1.77, |
||||
"Omega_deg": 131.78, |
||||
"omega_deg": 273.18, |
||||
"M_deg": -100.08, |
||||
}, |
||||
] |
||||
|
||||
AU = 1.496e11 # meters |
||||
|
||||
MOONS = [ |
||||
{ |
||||
"name": "Moon", |
||||
"mass": 7.35e22, |
||||
"radius": 1.738e6, |
||||
"parent": "Earth", |
||||
"a_km": 384400, |
||||
"e": 0.055, |
||||
"inc_deg": 5.16, |
||||
"Omega_deg": 125.08, |
||||
"omega_deg": 318.15, |
||||
"M_deg": 135.27, |
||||
}, |
||||
{ |
||||
"name": "Io", |
||||
"mass": 8.93e23, |
||||
"radius": 1.822e6, |
||||
"parent": "Jupiter", |
||||
"a_km": 421800, |
||||
"e": 0.004, |
||||
"inc_deg": 0.00, |
||||
"Omega_deg": 0.0, |
||||
"omega_deg": 49.1, |
||||
"M_deg": 330.9, |
||||
}, |
||||
{ |
||||
"name": "Europa", |
||||
"mass": 4.80e23, |
||||
"radius": 1.561e6, |
||||
"parent": "Jupiter", |
||||
"a_km": 671100, |
||||
"e": 0.009, |
||||
"inc_deg": 0.50, |
||||
"Omega_deg": 184.0, |
||||
"omega_deg": 45.0, |
||||
"M_deg": 345.4, |
||||
}, |
||||
{ |
||||
"name": "Ganymede", |
||||
"mass": 1.48e24, |
||||
"radius": 2.631e6, |
||||
"parent": "Jupiter", |
||||
"a_km": 1070400, |
||||
"e": 0.001, |
||||
"inc_deg": 0.20, |
||||
"Omega_deg": 58.5, |
||||
"omega_deg": 198.3, |
||||
"M_deg": 324.8, |
||||
}, |
||||
{ |
||||
"name": "Callisto", |
||||
"mass": 1.08e24, |
||||
"radius": 2.410e6, |
||||
"parent": "Jupiter", |
||||
"a_km": 1882700, |
||||
"e": 0.007, |
||||
"inc_deg": 0.30, |
||||
"Omega_deg": 309.1, |
||||
"omega_deg": 43.8, |
||||
"M_deg": 87.4, |
||||
}, |
||||
{ |
||||
"name": "Titan", |
||||
"mass": 1.35e24, |
||||
"radius": 2.575e6, |
||||
"parent": "Saturn", |
||||
"a_km": 1221900, |
||||
"e": 0.029, |
||||
"inc_deg": 0.30, |
||||
"Omega_deg": 78.6, |
||||
"omega_deg": 78.3, |
||||
"M_deg": 11.7, |
||||
}, |
||||
] |
||||
|
||||
|
||||
# Kepler conversion: M -> E -> nu |
||||
|
||||
def mean_to_true_anomaly(M_deg, e): |
||||
"""Convert mean anomaly (degrees) to true anomaly (radians) via Kepler's equation.""" |
||||
M = math.radians(M_deg) |
||||
E = solve_kepler_elliptical(M, e) |
||||
# tan(nu/2) = sqrt((1+e)/(1-e)) * tan(E/2) |
||||
tan_half_e = math.tan(E / 2.0) |
||||
tan_half_nu = math.sqrt((1.0 + e) / (1.0 - e)) * tan_half_e |
||||
nu = 2.0 * math.atan(tan_half_nu) |
||||
return normalize_angle(nu) |
||||
|
||||
|
||||
# Print TOML config |
||||
|
||||
def print_toml(): |
||||
print("# Moon Orbits Test Configuration") |
||||
print("# Auto-generated by scripts/precalc_moon_orbits.py") |
||||
print("# Data source: docs/planetary_data.md (JPL planetary facts)") |
||||
print("# Mean anomaly converted to true anomaly via Kepler's equation") |
||||
print() |
||||
|
||||
# Sun |
||||
print('[[bodies]]') |
||||
print('name = "Sun"') |
||||
print(f"mass = {SUN_MASS}") |
||||
print(f"radius = {SUN_RADIUS}") |
||||
print("parent_index = -1") |
||||
print('color = { r = 1.0, g = 1.0, b = 0.0 }') |
||||
print("orbit = { semi_major_axis = 0.0, eccentricity = 0.0, true_anomaly = 0.0 }") |
||||
print() |
||||
|
||||
# Planets |
||||
for p in PLANETS: |
||||
a_m = p["a_au"] * AU |
||||
inc = math.radians(p["inc_deg"]) |
||||
Omega = math.radians(p["Omega_deg"]) |
||||
omega = math.radians(p["omega_deg"]) |
||||
nu = mean_to_true_anomaly(p["M_deg"], p["e"]) |
||||
|
||||
print('[[bodies]]') |
||||
print(f'name = "{p["name"]}"') |
||||
print(f'mass = {p["mass"]}') |
||||
print(f'radius = {p["radius"]}') |
||||
print(f'parent_index = {p["parent"]}') |
||||
print('color = { r = 0.5, g = 0.5, b = 0.5 }') |
||||
print("orbit = {") |
||||
print(f" semi_major_axis = {a_m:.6e},") |
||||
print(f" eccentricity = {p['e']},") |
||||
print(f" inclination = {inc:.15f},") |
||||
print(f" longitude_of_ascending_node = {Omega:.15f},") |
||||
print(f" argument_of_periapsis = {omega:.15f},") |
||||
print(f" true_anomaly = {nu:.15f}") |
||||
print("}") |
||||
print() |
||||
|
||||
# Moons |
||||
for m in MOONS: |
||||
a_m = m["a_km"] * 1000.0 |
||||
inc = math.radians(m["inc_deg"]) |
||||
Omega = math.radians(m["Omega_deg"]) |
||||
omega = math.radians(m["omega_deg"]) |
||||
nu = mean_to_true_anomaly(m["M_deg"], m["e"]) |
||||
|
||||
print('[[bodies]]') |
||||
print(f'name = "{m["name"]}"') |
||||
print(f'mass = {m["mass"]}') |
||||
print(f'radius = {m["radius"]}') |
||||
parent_idx = {"Earth": 2, "Jupiter": 4, "Saturn": 5}[m["parent"]] |
||||
print(f'parent_index = {parent_idx}') |
||||
print('color = { r = 0.7, g = 0.7, b = 0.7 }') |
||||
print("orbit = {") |
||||
print(f" semi_major_axis = {a_m:.6e},") |
||||
print(f" eccentricity = {m['e']},") |
||||
print(f" inclination = {inc:.15f},") |
||||
print(f" longitude_of_ascending_node = {Omega:.15f},") |
||||
print(f" argument_of_periapsis = {omega:.15f},") |
||||
print(f" true_anomaly = {nu:.15f}") |
||||
print("}") |
||||
print() |
||||
|
||||
|
||||
# Print computed values summary (for verification) |
||||
|
||||
def print_summary(): |
||||
print("# === Computed True Anomalies ===") |
||||
print() |
||||
|
||||
for m in MOONS: |
||||
nu = mean_to_true_anomaly(m["M_deg"], m["e"]) |
||||
nu_deg = nu * 180.0 / math.pi |
||||
a_m = m["a_km"] * 1000.0 |
||||
mu = G * eval(f"{m['parent']}_MASS") if m["parent"] in globals() else 0 |
||||
|
||||
# Compute period |
||||
parent_mass = {"Earth": 5.97e24, "Jupiter": 1.898e27, "Saturn": 5.683e26}[m["parent"]] |
||||
mu = G * parent_mass |
||||
T = 2.0 * math.pi * math.sqrt(a_m**3 / mu) |
||||
T_days = T / 86400.0 |
||||
|
||||
print( |
||||
f'{m["name"]:10s}: M={m["M_deg"]:7.2f}deg -> nu={nu_deg:7.2f}deg ' |
||||
f"a={m['a_km']:>8.0f}km e={m['e']:.3f} " |
||||
f"T={T_days:.3f}d" |
||||
) |
||||
|
||||
print() |
||||
print("# === Initial positions (from true anomaly) ===") |
||||
print("# Format: name, r (m), nu (deg)") |
||||
|
||||
parent_masses = { |
||||
"Earth": 5.97e24, |
||||
"Jupiter": 1.898e27, |
||||
"Saturn": 5.683e26, |
||||
} |
||||
|
||||
for m in MOONS: |
||||
a_m = m["a_km"] * 1000.0 |
||||
nu = mean_to_true_anomaly(m["M_deg"], m["e"]) |
||||
pm = parent_masses[m["parent"]] |
||||
el = OrbitalElements( |
||||
a=a_m, e=m["e"], nu=nu, |
||||
inc=math.radians(m["inc_deg"]), |
||||
Omega=math.radians(m["Omega_deg"]), |
||||
omega=math.radians(m["omega_deg"]), |
||||
) |
||||
pos, vel = orbital_to_cartesian(el, pm) |
||||
r = vmag(pos) |
||||
print(f'{m["name"]:10s}: r={r:.3f} m, nu={nu*180/math.pi:.2f}deg') |
||||
|
||||
|
||||
if __name__ == "__main__": |
||||
print_summary() |
||||
print() |
||||
print("=" * 60) |
||||
print("# TOML CONFIG (copy below this line)") |
||||
print("=" * 60) |
||||
print() |
||||
print_toml() |
||||
Loading…
Reference in new issue