From fe6f9dcd5e6dd428cad74566714476ab5f12c0c5 Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Mon, 4 May 2026 12:24:51 -0400 Subject: [PATCH] remove extra decorative comment blocks in python scripts --- .../precalc_cartesian_to_elements_advanced.py | 16 - .../precalc_cartesian_to_elements_basic.py | 2 - scripts/precalc_extreme_eccentricity.py | 8 - scripts/precalc_extreme_timescales.py | 12 - scripts/precalc_inclined_orbits.py | 4 - scripts/precalc_moon_orbits.py | 337 ++++++++++++++++++ scripts/sim_engine.py | 28 -- 7 files changed, 337 insertions(+), 70 deletions(-) create mode 100644 scripts/precalc_moon_orbits.py diff --git a/scripts/precalc_cartesian_to_elements_advanced.py b/scripts/precalc_cartesian_to_elements_advanced.py index e839d76..80cf6d6 100644 --- a/scripts/precalc_cartesian_to_elements_advanced.py +++ b/scripts/precalc_cartesian_to_elements_advanced.py @@ -39,9 +39,7 @@ def report(name, original, recovered, fields): err = abs(orig_val - rec_val) print(f" {field:20s} = {orig_val:20.15e} -> {rec_val:20.15e} error = {err:.2e}") -# ============================================================================= # SECTION: eccentricity spectrum -# ============================================================================= print("=" * 70) print("SECTION: eccentricity spectrum: circular to highly hyperbolic") print("=" * 70) @@ -105,9 +103,7 @@ print(f"\n [8] Highly hyperbolic (e=10.0):") print(f" ecc error = {abs(rec_high_hyper.e - 10.0):.2e} (test tol: 1e-3)") print(f" a error = {abs(rec_high_hyper.a - (-1.0e10)):.2e} (test tol: 1e-2)") -# ============================================================================= # SECTION: inclination -# ============================================================================= print("\n" + "=" * 70) print("SECTION: inclination: zero, polar, and retrograde") print("=" * 70) @@ -133,9 +129,7 @@ rec_retro = roundtrip(retro) print(f"\n [3] Retrograde (inc=180 deg):") print(f" inc error = {abs(rec_retro.inc - math.pi):.2e} (test tol: 1e-4)") -# ============================================================================= # SECTION: true anomaly at key orbital positions -# ============================================================================= print("\n" + "=" * 70) print("SECTION: true anomaly at key orbital positions") print("=" * 70) @@ -158,9 +152,7 @@ for i, (nu_in, nu_exp, label) in enumerate(nu_tests): print(f" nu error = {nu_err:.2e} (test tol: 1e-6)") print(f" ecc error = {e_err:.2e} (test tol: 1e-4)") -# ============================================================================= # SECTION: quadrature at various eccentricities -# ============================================================================= print("\n" + "=" * 70) print("SECTION: quadrature at various eccentricities") print("=" * 70) @@ -177,9 +169,7 @@ for i, (e, e_tol, nu_tol) in enumerate(e_tests): print(f" a error = {a_err:.2e} (test tol: 1e-2)") print(f" nu error = {nu_err:.2e} (test tol: {nu_tol:.0e}) {'PASS' if nu_err <= nu_tol else 'FAIL'}") -# ============================================================================= # SECTION: large true anomaly values -# ============================================================================= print("\n" + "=" * 70) print("SECTION: large true anomaly values") print("=" * 70) @@ -201,9 +191,7 @@ for i, (nu_in, nu_exp, tol, label) in enumerate(large_nu_tests): print(f" a error = {a_err:.2e} (test tol: 1e-2)") print(f" nu error = {nu_err:.2e} (test tol: {tol:.0e}) {'PASS' if nu_err <= tol else 'FAIL'}") -# ============================================================================= # SECTION: 3D orientation with quadrature point -# ============================================================================= print("\n" + "=" * 70) print("SECTION: 3D orientation with quadrature point") print("=" * 70) @@ -217,9 +205,7 @@ print(f" inc error = {abs(rec.inc - math.pi/3.0):.2e} (test tol: 1e-4)") print(f" Omega error = {abs(rec.Omega - math.pi/4.0):.2e} (test tol: 1e-4)") print(f" omega error = {abs(rec.omega - math.pi/6.0):.2e} (test tol: 1e-4)") -# ============================================================================= # SECTION: multiple true anomaly points in sequence -# ============================================================================= print("\n" + "=" * 70) print("SECTION: multiple true anomaly points in sequence") print("=" * 70) @@ -236,9 +222,7 @@ for i, nu in enumerate(nu_seq): print(f" a error = {a_err:.2e} (test tol: 1e-2)") print(f" nu error = {nu_err:.2e} (test tol: 1e-6)") -# ============================================================================= # SECTION: hyperbolic orbit at quadrature point -# ============================================================================= print("\n" + "=" * 70) print("SECTION: hyperbolic orbit at quadrature point") print("=" * 70) diff --git a/scripts/precalc_cartesian_to_elements_basic.py b/scripts/precalc_cartesian_to_elements_basic.py index b821fd8..1cb1614 100644 --- a/scripts/precalc_cartesian_to_elements_basic.py +++ b/scripts/precalc_cartesian_to_elements_basic.py @@ -12,9 +12,7 @@ import sys, math sys.path.insert(0, 'scripts') from sim_engine import orbital_to_cartesian, cartesian_to_orbital_elements, vmag, OrbitalElements, G -# ============================================================================= # Test configuration: moderate eccentricity, zero inclination -# ============================================================================= mu = G * 5.972e24 a = 1.5e7 e = 0.5 diff --git a/scripts/precalc_extreme_eccentricity.py b/scripts/precalc_extreme_eccentricity.py index 1b5712a..04e8ca0 100644 --- a/scripts/precalc_extreme_eccentricity.py +++ b/scripts/precalc_extreme_eccentricity.py @@ -12,9 +12,7 @@ import sys, math sys.path.insert(0, 'scripts') from sim_engine import orbital_to_cartesian, cartesian_to_orbital_elements, vmag, OrbitalElements, G -# ============================================================================= # Spacecraft 0: Highly_Elliptical (e=0.99, a=6.5e8) -# ============================================================================= mu = G * 5.972e24 a0 = 6.5e8 e0 = 0.99 @@ -54,9 +52,7 @@ print(f"# v = {v0_pi:.6f} m/s") print(f"# dr = {abs(r0_pi - expected_r_apo0):.2e} m") print() -# ============================================================================= # Spacecraft 1: Near_Parabolic (e=0.99, a=7.0e8) -# ============================================================================= a1 = 7.0e8 e1 = 0.99 nu1 = 0.0 @@ -87,9 +83,7 @@ print(f"# dr_apo = {abs(r1_pi - expected_r_apo1):.2e} m") print(f"# v_peri > v_apo: {v1 > v1_pi}") print() -# ============================================================================= # Spacecraft 2: Slightly_Hyperbolic (e=1.05, a=-1.3e8) -# ============================================================================= a2 = -1.3e8 e2 = 1.05 nu2 = 0.0 @@ -116,9 +110,7 @@ print(f"# circular_vel = {circular_vel:.6f} m/s") print(f"# a < 0: {a2 < 0}") print() -# ============================================================================= # Velocity at different true anomalies for each spacecraft -# ============================================================================= print("# Velocity magnitudes at different true anomalies:") print("# (vis-viva: v = sqrt(mu * (2/r - 1/a)))") print() diff --git a/scripts/precalc_extreme_timescales.py b/scripts/precalc_extreme_timescales.py index e5f90d3..e158129 100644 --- a/scripts/precalc_extreme_timescales.py +++ b/scripts/precalc_extreme_timescales.py @@ -25,16 +25,12 @@ def circular_velocity(a, parent_mass): mu = G * parent_mass return math.sqrt(mu / a) -# =========================================================================== # Body definitions (from TOML) -# =========================================================================== earth_mass = 5.972e24 earth_radius = 6.371e6 sun_mass = 1.989e30 -# =========================================================================== # Spacecraft definitions and calculations -# =========================================================================== spacecraft = [ { "name": "Fast_Orbit_LEO", @@ -134,9 +130,7 @@ for sc in spacecraft: print() -# =========================================================================== # Geosynchronous period check -# =========================================================================== geo_a = 4.2164e7 geo_period = orbital_period(geo_a, earth_mass) sidereal_day_hours = 23.93447 @@ -149,9 +143,7 @@ print(f"# Period error: {abs(geo_period_hours - sidereal_day_hours):.6f} hours") print(f"# Period error: {abs(geo_period - sidereal_day_seconds):.6f} seconds") print() -# =========================================================================== # Jupiter-like 10-year propagation -# =========================================================================== jupiter_sc = spacecraft[2] jupiter_a = jupiter_sc["a"] jupiter_mu = G * jupiter_sc["parent_mass"] @@ -167,9 +159,7 @@ print(f"# Expected orbits = {expected_orbits:.6f}") print(f"# Expected true anomaly change = {expected_mean_anomaly % (2*math.pi):.10f} rad") print() -# =========================================================================== # Period consistency test: Mercury-like from different starting true anomalies -# =========================================================================== mercury_sc = spacecraft[1] mercury_a = mercury_sc["a"] mercury_e = mercury_sc["e"] @@ -187,9 +177,7 @@ for nu0_deg in [0, 90, 180, 270]: print(f"# After 1 period: true anomaly should return to {nu0_deg} deg") print() -# =========================================================================== # Low altitude orbit: check altitude above surface -# =========================================================================== low_sc = spacecraft[3] low_a = low_sc["a"] low_altitude = low_a - earth_radius diff --git a/scripts/precalc_inclined_orbits.py b/scripts/precalc_inclined_orbits.py index cc0b8a4..e50236d 100644 --- a/scripts/precalc_inclined_orbits.py +++ b/scripts/precalc_inclined_orbits.py @@ -13,9 +13,7 @@ import sys, math sys.path.insert(0, 'scripts') from sim_engine import orbital_to_cartesian, vmag, OrbitalElements, G -# ============================================================================= # Molniya orbit -# ============================================================================= a = 26540000.0 e = 0.74 inc = 1.107 @@ -40,9 +38,7 @@ print(f"#") print(f"# Period: {T:.6f} s = {T/3600:.6f} hours") print(f"# Half period: {T_half:.6f} s = {T_half/3600:.6f} hours") -# ============================================================================= # Generic inclined orbit -# ============================================================================= a2 = 10000000.0 e2 = 0.5 inc2 = math.radians(45) diff --git a/scripts/precalc_moon_orbits.py b/scripts/precalc_moon_orbits.py new file mode 100644 index 0000000..9eaf96e --- /dev/null +++ b/scripts/precalc_moon_orbits.py @@ -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() diff --git a/scripts/sim_engine.py b/scripts/sim_engine.py index 37eda77..2a9655a 100644 --- a/scripts/sim_engine.py +++ b/scripts/sim_engine.py @@ -17,9 +17,7 @@ from dataclasses import dataclass, field, replace from typing import Dict, Tuple, Any -# ============================================================================= # Constants -# ============================================================================= G = 6.67430e-11 PARABOLIC_TOLERANCE = 1e-3 @@ -28,9 +26,7 @@ KEPLER_MAX_ITER = 50 VEL_DRIFT_THRESHOLD = 1e-6 # m/s -# ============================================================================= # Vector operations -# ============================================================================= def vadd(a, b): return (a[0]+b[0], a[1]+b[1], a[2]+b[2]) @@ -68,9 +64,7 @@ def normalize_angle(angle): return angle -# ============================================================================= # Data structures -# ============================================================================= @dataclass class OrbitalElements: @@ -129,9 +123,7 @@ class Event: data: Dict[str, Any] = field(default_factory=dict) -# ============================================================================= # Burn direction vectors (local frame) -# ============================================================================ def get_burn_direction(direction, local_pos, local_vel): """Calculate burn direction vector in local frame.""" @@ -171,9 +163,7 @@ def apply_custom_burn(craft, delta_v_vec): craft.global_vel = vadd(craft.global_vel, delta_v_vec) -# ============================================================================= # Kepler equation solvers (exact C++ logic) -# ============================================================================= def get_initial_trial_value(mean_anomaly, eccentricity): """Initial guess for Kepler solver (C++ get_initial_trial_value).""" @@ -193,9 +183,7 @@ def solve_kepler_elliptical(mean_anomaly, eccentricity): return E -# ============================================================================= # Coordinate transforms -# ============================================================================= def orbital_to_cartesian(elements, parent_mass): """Convert orbital elements to local position/velocity vectors.""" @@ -352,9 +340,7 @@ def cartesian_to_orbital_elements(pos, vel, parent_mass): return elements -# ============================================================================= # Propagation -# ============================================================================= def propagate(elements, dt, parent_mass): """Propagate orbital elements forward by dt. Returns new elements.""" @@ -399,9 +385,7 @@ def propagate(elements, dt, parent_mass): raise NotImplementedError("hyperbolic propagation not yet implemented") -# ============================================================================= # Global coordinate computation -# ============================================================================= def compute_global_coordinates(bodies): """ @@ -418,9 +402,7 @@ def compute_global_coordinates(bodies): body.global_vel = vadd(body.local_vel, parent.global_vel) -# ============================================================================= # Velocity drift check -# ============================================================================= def check_velocity_drift(body, parent, parent_mass): """ @@ -437,9 +419,7 @@ def check_velocity_drift(body, parent, parent_mass): body.orbit = cartesian_to_orbital_elements(body.local_pos, body.local_vel, parent_mass) -# ============================================================================= # Body physics update -# ============================================================================= def update_body(bodies, body_index, dt): """ @@ -458,9 +438,7 @@ def update_body(bodies, body_index, dt): body.local_pos, body.local_vel = orbital_to_cartesian(body.orbit, parent.mass) -# ============================================================================= # Spacecraft physics update -# ============================================================================ def update_spacecraft(spacecraft_list, bodies, dt): """ @@ -492,9 +470,7 @@ def compute_global_coordinates_spacecraft(spacecraft_list, bodies): craft.global_vel = vadd(craft.local_vel, parent.global_vel) -# ============================================================================= # TOML config loader -# ============================================================================= def load_config(config_path): """Load a TOML 1.0 config file and return parsed data.""" @@ -605,9 +581,7 @@ def initialize_spacecraft(spacecraft_list, bodies): craft.global_vel = (0.0, 0.0, 0.0) -# ============================================================================= # Initialization -# ============================================================================= def initialize_bodies(bodies): """ @@ -629,9 +603,7 @@ def initialize_bodies(bodies): body.global_vel = (0.0, 0.0, 0.0) -# ============================================================================= # Simulator — public API -# ============================================================================= class Simulator: """