For e < 1e-10, compute argument of latitude instead of hardcoding
true_anomaly = 0.0. This preserves position consistency after circularization
burns.
- Move ascending node vector computation earlier so it's available for
the circular orbit case
- For well-defined ascending nodes (sin_i > 1e-6): compute argument of
latitude using the ascending node direction and h×n perpendicular
- For nearly-coplanar orbits (sin_i <= 1e-6): use atan2(y, x) directly
to avoid numerically unstable ascending node direction
For coplanar orbits (inclination < 0.01 rad), compute omega from the
eccentricity vector (longitude of periapsis) instead of forcing omega=0.
This is necessary for correct post-burn orbital element conversion,
critical for Hohmann transfer accuracy.
Update test_omega_debug to accept the new behavior: omega is computed
from the eccentricity vector direction for coplanar orbits.
- Fix omega = π bug for near-zero inclination orbits (src/orbital_mechanics.cpp:277)
- Added inclination threshold (0.01 rad) before using atan2 path
- Prevents unstable omega calculation when h_vec.y is tiny
- Fix true anomaly trigger firing at wrong location (src/maneuver.cpp:172-213)
- Replaced binary search with O(1) analytical time calculation
- Uses mean motion: n = √(μ/a³) to compute exact time to target
- Added true_anomaly_to_eccentric_anomaly() function
See docs/planning/periapsis_burn_bug_analysis.md for technical details
See docs/periaapsis_burn_test_results.md for test results and next steps
Test status:
- Issue 1 (omega): ✅ FIXED and validated
- Issue 2 (triggers): ✅ IMPLEMENTED, 4 tests fail due to design (expected behavior)
- Full suite: 139 passed, 4 failed (all periapsis burn tests)
Next steps: Decide on test tolerance approach (see docs/periaapsis_burn_test_results.md)
Added planning document detailing two separate issues causing
periapsis burns to execute at incorrect locations:
Issue 1: Omega = π instead of 0 for near-zero inclination orbits
- Unstable atan2 calculation in cartesian_to_orbital_elements()
- Occurs when n_mag > 1e-10 due to tiny h_vec.y from i ≈ 0
- Fix: Add inclination threshold (0.01 rad) before using atan2 path
Issue 2: True anomaly trigger fires early at wrong location
- angle_between() detects upcoming crossing, fires immediately
- Executes burn at current position instead of waiting for target angle
- Proposed fixes: Defer execution, interpolate, or remove future check
Also cleaned up:
- Removed temporary debug files (debug_test_burn.cpp, debug_trace.cpp)
- Removed debug output from src/orbital_mechanics.cpp
- Kept tests/test_omega_debug.cpp as it validates Issue 1 fix
See docs/planning/periapsis_burn_bug_analysis.md for full details.
- Function was never called anywhere in the codebase
- Functionality now handled by propagate_orbital_elements() with proper Newton-Raphson iteration
- Removes dead code and simplifies the API
- Switch spacecraft and bodies from RK4 to analytical propagation using propagate_orbital_elements()
- Fixed three bugs in hyperbolic orbit propagation (formula errors and insufficient Newton-Raphson iterations)
- Add orbital element reconstruction after burns and SOI transitions
- Zero energy drift for circular, elliptical, parabolic, and hyperbolic orbits
- All 132 tests passing (240,294 assertions)
- Add solve_barker_equation() function using cubic formula: D + D³/3 = M
- Integrate Barker's equation into propagate_orbital_elements() for parabolic orbits
- Add comprehensive test suite (11 tests, 239 assertions) following TDD
- Use cbrt() for cube root (handles negative numbers properly)
- Parabolic propagation now uses exact analytical solution instead of iterative solver
- All 93 tests passing (239,872 assertions)
- Add 16 minimal comments documenting formulas in orbital_mechanics.cpp
- Fix parabolic test to use orbital_elements_to_cartesian() instead of manual velocity
- Tighten parabolic test tolerance from 1e10 to 1e3 (7 orders of magnitude)
- Reduce parabolic test error from 6.5% to machine precision
- Add near-parabolic detection for |1-e| < 0.01
- Use cosine/sine formulation with atan2() to preserve quadrant
- Avoid catastrophic cancellation in sqrt((1+e)/(1-e)) factor
- Clamp values to [-1,1] to handle precision issues
- Maintain original formula for well-behaved orbits
- Add PARABOLIC_TOLERANCE = 1e-3 constant for consistent detection
- Replace inconsistent thresholds (0.005, 0.98, 1.02) across 5 files
- Refactor orbital_elements_to_cartesian() to use semi-latus rectum as primary parameter
- Eliminate 3 separate code branches with unified formulas for all orbit types
- Improve numerical stability for parabolic and near-parabolic orbits
- Reduce code complexity: -23 lines net
- Convert 64 test assertions from Approx() to WithinAbs() in 2 new test files
- Add WithinAbs() testing guidelines to AGENTS.md
- Fix cartesian_to_orbital_elements(): eccentricity vector calculation,
true anomaly normalization, parabolic semi-latus rectum handling
- Add 2 new test files for edge cases and quadrature points
New modular API:
- solve_kepler_elliptical(M, e): Newton-Raphson for E - e·sin(E) = M
- solve_kepler_hyperbolic(M, e): Solver for H - e·sinh(H) = M
- eccentric_to_true_anomaly(E, e): Convert eccentric to true anomaly
- hyperbolic_to_true_anomaly(H, e): Convert hyperbolic to true anomaly
- mean_anomaly_to_true_anomaly(M, e): Unified wrapper (dispatches based on e)
Changes:
- Renamed solve_kepler_equation() → solve_kepler_elliptical() for clarity
- Extracted KEPLER_TOLERANCE and KEPLER_MAX_ITERATIONS constants
- Separated hyperbolic solver logic from combined function
- Fixed test_newton_raphson_convergence to verify Kepler's equation
(instead of incorrectly expecting E ≈ M for small e)
- Added TODO comment for future cartesian_to_orbital_elements refactoring
Bug: r_dot_e was incorrectly divided by mu instead of r_mag*e_mag,
causing cos_nu values of ~7.5e6 (far outside [-1,1]) which made acos() return NaN.
Fix: Calculate cos_nu correctly as r_dot_e/(r_mag*e_mag) and clamp
to [-1, 1] to handle floating-point precision errors.
This fixes all cartesian <-> orbital elements round-trip conversion tests.
Added 6 test files for Newton-Raphson solver and analytical propagation:
- test_cartesian_to_elements_basic.cpp: Tests state vector ↔ orbital elements conversion
- test_newton_raphson_convergence.cpp: Tests Newton-Raphson solver convergence behavior
- test_analytical_propagation_apsides.cpp: Tests propagation through orbital apsides
- test_analytical_propagation_timesteps.cpp: Tests propagation with various timesteps
- test_extreme_eccentricity.cpp: Tests near-parabolic and hyperbolic orbits
- test_precision_boundaries.cpp: Tests exact boundary value handling
Implemented core orbital mechanics functions:
- solve_kepler_equation(): Newton-Raphson solver with 1e-10 tolerance
- get_initial_trial_value(): Series expansion initial guess
- cartesian_to_orbital_elements(): State vectors to orbital elements conversion
- propagate_orbital_elements(): Analytical propagation using Kepler's equation
Updated test plan document with current progress and remaining tests.
Test status: 66 passed, 14 failed (out of 80 test cases)
- Failing tests are expected: implementation needs debugging
- Config validation issues fixed by adjusting orbital parameters
- Modify orbital_elements_to_cartesian() to apply z-x-z Euler rotations
- Apply rotation: R_z(Ω) · R_x(i) · R_z(ω) to position and velocity
- Fix generic inclined orbit test: set argument_of_periapsis to π/2
- Molniya position tests now pass with non-zero z-coordinates
- Period test still fails (orbit tracker doesn't handle 3D orbits yet)
- Rewrote orbital_mechanics.cpp with complete velocity component handling
- Added proper parabolic case in velocity section using semi_latus_rectum
- Parabolic velocity: vx = -sqrt(μ/p)*sin(ν), vy = sqrt(μ/p)*(1+cos(ν))
- Updated parabolic_comet.toml to use semi_latus_rectum = 2.992e11 (p=2 AU)
- With true_anomaly = 0.0, comet at perihelion (r = 1 AU, v = 42.13 km/s)
Test Results:
- 31 passed / 1 failed (up from 29/3)
- Parabolic orbit tests: All 3 assertions pass ✓
- Position and velocity mathematically correct for parabolic orbits
- Total energy approximately zero (FP precision: ~200 J/kg)
Note: test_invalid_parent_assignment still fails - pre-existing issue unrelated to parabolic implementation.
- Added union in OrbitalElements (semi_major_axis | semi_latus_rectum)
- Updated config_loader to parse and validate semi_latus_rectum for parabolic orbits
- Validation ensures correct parameter (semi_latus_rectum for e≈1, semi_major_axis for other)
- Updated orbital_mechanics.cpp parabolic case: r = p / (1 + cos(ν))
- Updated parabolic_comet.toml to use semi_latus_rectum = 1.496e11
- Tolerance for parabolic detection: |e - 1.0| < 0.005
Note: Parabolic test still fails velocity tolerance check due to floating-point precision
at r=0.5 AU, E ≈ -200 J/kg (not exactly 0)
- Added orbital_mechanics module with orbital_elements_to_cartesian()
- Updated initialize_orbital_objects() to use orbital mechanics
- Added validate_initial_positions() for post-initialization checking
- Fixed test files to use global_position/global_velocity
- Updated config loader to support spacecraft altitude parameter
- Fixed orbital_mechanics.cpp velocity calculation bug (removed duplicate scaling)
- Updated Makefile to include orbital_mechanics.o in test build
- Renamed simulation.h OrbitalMetrics to OrbitalAnalysis to avoid conflict
- Added docs/parabolic_union_implementation.md for parabolic orbit support plan
Note: Test configs still need manual fix for orbit table TOML syntax
cinnaboot